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1.  INTRODUCTION 


It  has  been  observed  experimentally  that  jet  fuels,  when  heated,  undergo  chemical  reactions  that 
eventually  result  in  sediment/ deposit  formation.  In  an  aircraft  engine,  under  normal  operating  condi¬ 
tions,  the  temperature  in  the  fuel  supply  systems,  such  as  the  heat  exchanger,  fuel  lines,  and  injection 
nozzles,  is  typically  in  the  range  of  100-300*C.  At  these  temperatures,  the  fuel  tends  to  undergo 
degradation  reactions,  resulting  in  the  formation  of  deposits.  This  can  adversely  affect  the  operation 
of  the  aircraft  by  clogging  the  fuel  lines  and  injection  nozzles.  With  the  increased  use  of  advanced 
turbines  designed  to  operate  at  higher  temperatures,  thermal  degradation  problems  are  expected  to 
increase.  Also,  in  supersonic  aircraft,  a  considerable  increase  in  the  metal  skin  temperature  is  possible 
due  to  aerodynamic  heating.  For  example,  it  has  been  estimated  that  for  an  aircraft  cruising  at  Mach 
2.7,  the  external  skin  temperature  can  reach  into  the  range  of  450-500‘F  (230-260*C),  and  the  tem¬ 
perature  of  an  uninsulated  fuel  tank  could  rise  to  430*F  (220*C).  Hydrocarbon  fuel  exposed  to  such 
temperature  stress  can  undergo  degradation  and  form  deposits  (Taylor,  1969a).  Our  knowledge  of  the 
chemical  structure  of  fuel  deposits,  the  mechanisms  for  deposit  formation,  and  the  locations  where  the 
deposits  collect  is  currently  limited.  By  understanding  the  processes  involved,  suitable  modifications 
to  the  fuel,  handling  techniques,  or  surface  material  can  be  made  to  eliminate  or  minimize  deposit  for¬ 
mation.  This  would  result  in  improvements  in  the  design  of  fuel  systems  for  future  aircraft. 

The  susceptibility  of  hydrocarbon  fuel  to  thermal  degradation  is  usually  called  thermal  instability 
and  is  influenced  by  many  factors,  such  as  fuel  type,  temperature,  pressure,  and  the  material  in  con¬ 
tact  with  the  fuel  (Marteney,  1988;  Marteney  and  Spadaccini,  1986;  Taylor,  1967,  1969b;  CRC,  1978, 
1979;  Dahlin  et  al.,  1981;  Daniel,  1983,  1985).  The  effects  of  these  factors  on  the  formation  of  depo¬ 
sits  have  been  investigated  experimentally  under  both  nonflowing  and  flowing  conditions.  In  many  of 
these  investigations  (e.g.,  Cohen,  1980;  Taylor,  1969b,  1970;  Taylor  and  Frankenfeld,  1978),  only  the 
overall  rate  of  deposit  was  measured.  The  chemical  structure  of  fuel  deposits  has  been  determined  in 
some  experiments  for  various  types  of  jet  fuels  (e.g.,  Hemlick  and  Seng,  1984;  Nixon  and  Henderson, 
1966).  The  presence  of  molecular  oxygen  in  fuel  lines  appears  to  play  a  major  role  in  the  chemical 
degradation  of  the  fuel  (e.g.,  Nixon  and  Henderson,  1966;  Taylor,  1974)  and  recent  experiments  indi¬ 
cate  that  autooxidation  reactions  are  the  primary  process  of  deposit  formation  at  temperatures  less 
than  535K  (Marteney,  1988;  Marteney  and  Spadaccini,  1986).  However,  it  has  also  been  observed 
that  deoxygenated  fuel  forms  deposits  if  the  temperature  is  high  enough  (Taylor  and  Frankenfeld, 
1978).  Deposits  formed  under  these  conditions  differ  in  composition  from  those  formed  in  the  pres¬ 
ence  of  oxygen.  Furthermore,  in  deoxygenated  fuels,  the  addition  of  trace  compounds,  such  as  nitro¬ 
gen  compounds,  decreases  the  thermal  stability  (Nixon,  1962).  A  similar  decrease  in  stability  was 
observed  when  trace  amounts  of  sulphur  compounds,  i.e.,  thiols,  sulfides,  and  disulfides,  were  added 
to  deoxygenated  hydrocarbon  fuels  (Taylor,  1976;  Taylor  and  Wallace,  1968).  This  decrease  in  stabil¬ 
ity  was  attributed  to  a  process  of  free-radical-initiated  chain  reactions,  although  details  of  the  reaction 
mechanism  are  unknown.  Metal  surfaces  and  dissolved  metal  salts  can  also  decrease  thermal  stability 
(e.g.,  Schenk  et  al.,  1971;  Taylor,  1969b;  Kendall  and  Mills,  1986).  On  the  other  hand,  removal  of 
molecular  oxygen  from  the  fuel  has  shown  an  increase  in  stability,  resulting  in  reduced  deposit  forma¬ 
tion  (Taylor,  1974).  Thus,  the  process  of  thermal  degradation  of  both  fuels  with  molecular  oxygen  and 
deoxygenated  fuels  appears  to  be  influenced  by  many  factors,  some  of  which  are  still  unknown. 

The  formation  of  deposits  under  conditions  in  which  the  fuel  is  flowing  over  heated  surfaces  has 
also  been  extensively  studied.  Earlier  experiments  employed  the  small-scale  thermal  stability  test  dev¬ 
ices  such  as  the  Jet  Fuel  Thermal  Oxidation  Tester  or  JFTOT  (CRC,  1978;  UTC,  1988).  More 
recently,  an  extensive  test  program  was  conducted  by  Marteney  (1988)  and  Marteney  and  Spadaccini 
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(1986)  to  determine  the  thermal  stability  and  heat  transfer  characteristics  of  JP-5  and  several 
kerosene-type  fuels  using  small,  resistively  heated  tubes.  Data  under  various  operating  conditions  was 
obtained  with  a  high  degree  of  repeatability.  It  was  determined  that  the  surface  temperature  was  the 
key  parameter  affecting  the  deposit  formation,  although  additional  parameters  such  as  the  fuel  flow 
rate  and  test  duration  were  also  identified  as  important  physical  parameters.  One  of  the  important 
results  of  these  studies  was  that  the  thermal  degradation  process  occurs  in  two  stages,  each  with  its 
own  unique  activation  energy.  The  first  stage  occurs  at  relatively  low  temperatures  (less  than  535K) 
and  appears  to  be  primarily  autooxidation-limited,  whereas  the  second  stage  occurs  at  higher  tempera¬ 
tures  (535-61 5K)  and  appears  to  be  combination  of  autooxidation  and  hydrocarbon  pyrolysis.  How¬ 
ever,  in  another  set  of  experiments  (Chin  et  al.  1989),  with  heated  flowing  fuels  over  heated  metal  sur¬ 
faces,  the  deposition  characteristics  were  in  many  significant  respects  quite  different  than  Marteney’s 
data.  Thus,  there  are  still  some  unanswered  questions  regarding  the  exact  thermochemical  processes 
involved  in  the  thermal  degradation  of  jet  fuels. 

Since  experimental  studies  are  at  present  incapable  of  providing  a  detailed  understanding  of  the 
various  phenomena  occurring  during  thermal  degradation  of  jet  fuels,  there  is  a  need  to  develop  an 
alternative  approach  that  may  provide  a  fundamental  improvement  in  our  understanding  of  the  phy¬ 
sics  of  the  degradation  process.  One  such  alternative  approach  is  the  use  of  computational  fluid 
dynamics  (CFD)  to  develop  models  that  can  be  systematically  used  to  study  the  fluid  dynamics,  heat 
transfer,  and  chemical  processes  that  occur  during  thermal  degradation  of  jet  fuels.  Since  the  various 
parameters  can  be  varied  in  a  CFD  approach  relatively  easily  as  compared  to  the  experimental 
approach,  CFD  models  may  provide  the  necessary  flexibility  to  predict  the  sensitivity  of  the  parame¬ 
ters  that  govern  the  deposition  process.  However,  it  is  clear  that  to  properly  use  CFD  models,  they 
must  first  be  validated  by  comparing  the  predictions  with  experimental  data.  Therefore,  the  availabil¬ 
ity  of  accurate  experimental  data  is  essential  for  the  successful  development  of  CFD  models. 

Recently,  CFD  models  have  been  developed  by  Oh  et  al.  (1989),  Roquemore  et  al.  (1989),  and 
Krazinski  et  al.  (1990).  Three-dimensional  effects  were  investigated  by  Oh  et  al.  (1989),  who  found 
that  a  complex  spiraling  3D  flow  exists  in  the  JFTOT.  In  the  study  by  Krazinski  et  al.  (1990),  a  ther¬ 
mal  deposition  model  was  incorporated  into  an  existing  CFD  code  which  solved  the  Reynolds- 
averaged  equations  of  mass,  momentum,  and  energy.  Steady-state  solutions  were  obtained  using  an 
efficient  multigrid  technique.  The  thermal  deposition  process  was  modeled  by  three  global  Arrhenius 
reactions  in  which  the  thermal  decomposition  of  the  fuel  was  assumed  to  be  primarily  due  to  the 
autooxidation  reaction  with  the  dissolved  oxygen  in  the  jet  fuel.  Since  various  parameters  had  to  be 
calibrated,  the  experimental  data  of  Marteney  (1988)  and  Marteney  and  Spadaccini  (1986)  was  first 
employed  to  calibrate  the  unknown  parameters.  The  model  was  then  tested  on  other  fuel  decomposi¬ 
tion  data.  Good  agreement  was  obtained  in  these  studies  for  both  laminar  and  turbulent  fuel  flows. 

The  CFD  model  described  above,  however,  had  some  limitations,  some  of  which  are  currently 
being  addressed  (Roquemore,  private  communication).  For  example,  the  steady-state  analysis  did  not 
address  the  nonlinear  time-dependent  behavior  of  the  deposition  process.  Experiments  by  Marteney 
(1988)  showed  that  as  the  test  duration  was  increased,  the  deposition  rate  also  increased.  Such  a  non¬ 
linear,  time-dependent  phenomena  must  be  modeled  using  a  time-accurate  model.  Furthermore,  the 
model  did  not  account  for  the  formation  of  solid  particles  in  the  flowing  fuel  and  its  subsequent  tran¬ 
sport  and  impaction  on  the  wall.  As  the  deposit  grows  on  the  wall,  it  will  affect  the  heat  transfer 
between  the  heated  rod  and  the  fuel  (Marteney,  1988),  a  process  that  was  also  not  modeled  in  the 
above-mentioned  study.  In  spite  of  these  limitations,  the  results  obtained  so  far  are  quite  impressive 
in  demonstrating  the  capability  of  the  CFD  approach  to  model  and  predict  the  fuel  decomposition  and 
deposit  formation  process. 
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In  the  Phase  I  study  described  here,  a  mathematical  model  was  developed  that  takes  into  account 
the  time-dependent  interactions  among  the  heat  transfer  in  the  heating  element,  the  heat  transfer 
across  the  thin  layer  of  deposit,  the  heat  loss  to  the  flowing  fluid  and  the  chemical  kinetics  of  the  fuel 
degradation  reactions.  The  complex  coupling  among  the  unit  processes  is  modeled  by  a  system  of  dif¬ 
ferential  equations.  In  many  respects,  the  model  developed  here  is  similar  to  the  CFD  model 
developed  by  Roquemore  et  al.  (1989)  and  Krazinski  et  al.  (1989).  There  are  some  differences,  how¬ 
ever,  which  will  become  apparent  in  the  subsequent  sections.  For  demonstration  purposes,  a  simple 
fuel  flow  system  is  chosen  which  approximately  models  a  small-scale  test  device  (i.e.,  JFTOT)  that  has 
been  used  extensively  in  many  experimental  studies.  The  results  of  this  study  are  documented  in  this 
report. 


2.  TECHNICAL  OBJECTIVES 

The  objectives  of  the  Phase  I  research  were  as  follows: 

1.  To  construct  a  mathematical  model  of  fuel  deposition  in  a  simple  flowing  system  by  using 
existing  knowledge  on  each  unit  process. 

2.  To  demonstrate  the  mathematical  model  by  applying  it  to  a  simple  fuel  system  with  known 
experimental  results. 

3.  To  compare  the  results  with  experimental  measurements. 

As  described  in  this  report,  all  these  objectives  have  been  successfully  achieved. 


3.  FORMULATION  OF  THE  MODEL 

In  this  section,  we  describe  in  some  detail  the  formulation  of  the  model.  The  formulation  was  car¬ 
ried  out  in  a  general  sense  so  that  eventual  application  to  realistic  fuel  systems  could  be  carried  out  in 
the  next  phase.  However,  since  there  are  various  aspects  of  the  coupling  and  the  chemical  kinetics 
that  are  still  poorly  understood,  some  simplifications  to  the  full  model  were  made  before  carrying  out 
the  computations.  These  simplifications  will  be  discussed  later. 

Due  to  the  limited  resources  available  for  the  Phase  I  research,  we  investigated  a  fuel  flow  system 
with  a  simple  model.  The  fuel  system  modeled  is  the  JFTOT  test  device  configuration  used  in  the 
experiments.  It  is  well-known  that  the  flow  field  in  the  JFTOT  is  three-dimensional  (CRC,  1978, 
1979;  UTC,  1988;  Oh  et  al.,  1989);  however,  for  this  phase  of  study  we  restrict  ourselves  to  axisym- 
metric  geometry.  Extension  to  three-dimensional  flow  can  be  accomplished  in  the  next  phase.  The 
flow  configuration  considered  here  is  shown  in  Figure  1.  It  consists  of  a  circular  pipe  of  length  /  and 
internal  diameter  of  2r3,  with  a  heating  rod  of  diameter  2 rx  at  the  center.  The  rod  is  heated  by  an 
electrical  element  at  its  center  so  that  the  heating  process  can  be  simulated  by  a  prescribed  heat  flux 
condition  at  the  centerline.  This  is  similar  to  the  experimental  condition  of  Marteney  and  was  also 
used  in  the  CFD  model  developed  by  Roquemore  et  al.  (1989)  and  Krazinski  et  al.  (1989).  However, 
the  CFD  code  was  written  using  a  general  format  so  that  various  types  of  boundary  conditions  can  be 
implemented  easily  without  any  major  modifications. 
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Figure  1.  An  Axisymmetric  Fuel  Flow  System  with  Characteristic  JFTOT  Dimensions 
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There  are  five  main  unit  processes  that  must  be  modeled  and  coupled  together  to  simulate  the 
problem  of  deposit  formation  due  to  thermal  degradation  of  the  flowing  fuel: 

1.  Heat  transfer  in  the  heating  rod. 

2.  Heat  transfer  across  the  fuel  deposit. 

3.  Mass,  momentum,  and  heat  transfer  in  the  fluid  phase. 

4.  Solid  and  fuel-soluble  deposit  production  and  transport  in  the  fluid  phase. 

5.  Formation  of  the  solid  deposit  on  the  heated  wall  due  to  fuel  degradation  reactions. 

These  five  unit  processes  are  described  separately  in  the  following  subsections.  Note  that  the  third 
and  fourth  unit  processes  are  coupled  phenomena  in  the  fluid  phase  that  must  be  modeled  together. 


Nondimensionalization 

To  implement  the  model  into  a  computational  code,  it  is  preferable  that  the  governing  equations  be 
nondimensional.  Therefore,  all  the  equations  were  nondimensionalized  by  choosing  reference  quanti¬ 
ties.  We  chose  a  reference  length  L  which  for  the  present  study  was  the  radius  of  the  heating  rod  (rj. 
A  reference  velocity  u,  was  chosen  for  the  flowing  fuel  and  taken  to  be  the  mean  free  stream  velocity 
of  the  fuel  at  the  inlet.  The  reference  time  scale  is  thus  tT  =  L/ur.  The  reference  temperature  Tr  was 
chosen  to  be  the  temperature  of  the  fuel  at  the  inlet  (typically  room  temperature).  All  transport  pro¬ 
perties,  such  the  kinematic  viscosity  of  the  fuel  and  the  thermal  conductivities  of  the  rod  (/c,),  the 
deposit  (/cd),  and  the  fuel  ( Kf ),  are  nondimensionalized  by  their  respective  values  at  the  reference  tem¬ 
perature  (Tr). 

In  the  following  formulation,  the  equations  are  presented  primarily  in  the  dimensional  form, 
although  in  most  cases  the  effect  of  nondimensionalization  is  also  described.  All  dimensional  quanti¬ 
ties  have  units  consistent  with  the  metric  system. 


3.1  Heat  Transfer  in  the  Heating  Rod 

The  first  unit  process  that  needs  to  be  modeled  is  the  heat  transfer  process  in  the  heating  rod.  The 
model  for  the  heating  rod  consists  of  the  heat  source  at  a  given  location,  the  axial  and  radial  heat 
transfer  in  the  rod,  and  the  heat  loss  to  the  deposit  layer  due  to  conduction.  The  equation  governing 
the  temperature  distribution  in  the  rod  can  be  written  as 

dT, 

p,Cp,~aT  =  V(*' VFl)  (1) 

where  Tx  is  the  temperature  in  the  rod  and  p,,  and  Cp<  are  the  thermal  conductivity  coefficient, 
the  density  of  the  rod  material,  and  the  specific  heat  of  the  rod  material,  respectively.  In  the  JFTOT, 
the  rod  ij  made  of  stainless  steel.  For  the  present  study,  the  density  (pt)  and  the  specific  heat  (C,)  are 
assumed  to  be  constants  in  the  temperature  range  studied,  while  the  thermal  conductivity  of  the  stain¬ 
less  steel  rod  is  a  weak  function  of  the  temperature.  Thus,  k,  is  prescribed  by  a  second-order  polyno¬ 
mial  in  temperature  (i.e.,  k,  =  a  +  bT  +  cT 2),  where  the  coefficients  are  determined  by  fitting  the 
polynomial  to  the  data  in  the  temperature  range  of  100-500*C. 
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In  Equation  (1),  V  is,  in  general,  the  three-dimensional  differential  operator.  For  the  present 
study,  we  reduce  it  to  the  axisymmetric  differential  operator  so  that  after  nondimensionalizing,  the 
final  equation  for  transport  of  heat  in  the  rod  can  be  written  as 


-vr-ar 

a 


d 

_  dT 

K, - 

+  — — 

—  dT 

TK, - 

d X 

dx 

r  dr 

dr  J 

(2) 


where  the  bar  denotes  the  nondimensional  quantities  and,  for  simplicity,  the  subscript  I  in  Equation 
(1)  has  been  neglected.  Also,  a,  is  a  constant  that  appears  due  to  the  nondimensionalization  and  is 
defined  as 


o. 


Pf,  Cpfl  «», 

P;  Cp.t  Kfr 


1 

ReLPr, 


(3) 


where  the  rod  transport  properties  have  been  rescaled  using  the  reference  values  for  the  fuel.  Here, 
ReL  =  pjUrL/^ir  is  the  Reynolds  number  of  the  fuel  flow  based  on  the  reference  length  L  and  fuel 
properties  at  the  reference  temperature,  and  PrT  =  is  the  Prandtl  number  of  the  fuel  at  the 

reference  temperature.  The  physical  domain  in  which  this  equation  is  valid  is  defined  by  0  <  x  <  l  in 
the  axial  direction,  where  l  is  the  length  of  the  rod,  and  by  0  <  r  <  rx  in  the  radial  direction,  so  that 
when  nondimensionalized,  the  axial  domain  is  0  <  x  <  l/rx  and  the  radial  domain  is  0  <  r  <  1. 

Equation  (2)  is  solved  subject  to  proper  boundary  conditions.  The  formulation  allows  us  to  imple¬ 
ment  any  type  of  boundary  conditions  at  the  boundaries.  Therefore,  in  general,  the  dimensional  boun¬ 
dary  conditions  for  the  heating  rod  are 


dT , 

dx 

=  <7o 

at 

x  =  0 

(4) 

dTy 

dx 

=  Qi 

at 

x  =  / 

(5) 

dTx 

dr 

=  Qi 

at 

r  =  0 

(6) 

dT  x 

dr 

ll 

at 

r  =  r  i 

(7) 

Here,  the  subscript  1  denotes  the  variables  in  the  heating  rod,  and  the  subscript  2  denotes  those  in  the 
deposit  layer.  The  boundary  conditions  at  x  =  0  and  x  =  /  allows  the  possibility  of  heat  loss  (if 
q0  <  0  and  qi  <  0)  or  heat  addition  (if  q0  >  0  and  qt  >  0)  or  insulated  wall  (if  q0  =  qt  —  0).  There  is 
also  an  option  to  use  constant  temperature  end  conditions  [i.e.,  T(x  =  0)=  Tr  and  T(x  =  1)  =  Tr).  If 
the  rod  is  heated  by  a  heating  element  as  is  modeled  in  the  present  JFTOT  configuration,  then  at 
r  =  0,  the  heat  flux  is  prescribed  by  qt.  Finally,  the  boundary  condition  given  by  Equation  (7)  deter¬ 
mines  the  heat  loss  from  the  heating  rod  to  the  deposit  layer. 


TR-502/03-90 


6 


3.2  Heat  Transfer  Across  the  Fuel  Deposit 


The  second  unit  process  is  the  heat  transfer  across  the  deposit  layer.  The  deposit  is  defined  by  the 
axial  domain  0  <  x  <  l  and  the  radial  domain  of  rj  <  r  <  r2.  Note  that  the  deposit  thickness  will 
grow  in  the  axial  direction  and  in  time  and  thus,  r2  =  r2(x,  /).  Since  the  deposit  layer  is  usually  very 
thin,  the  heat  transfer  in  the  deposit  layer  can  be  modeled  as  a  steady-state  process,  i.e.,  the  equation 
governing  the  heat  conduction  through  the  layer  is  then 

V-(k  dvr2)  =  0  (8) 


Furthermore,  since  the  layer  is  very  thin,  the  axial  variation  of  temperature  in  the  deposit  can  be 
neglected  when  compared  to  its  radial  variation.  Thus,  Equation  (8)  reduces  to  a  one-dimensional 
form: 


d_ 

dr 


=  0 


(9) 


This  equation  can  be  directly  integrated  to  obtain 


T?  *5  Q 

/  KjdT  =  /— dr  (10) 

7l+  f!  r 


where  Cd  is  the  as  yet  unknown  constant  of  integration,  Tx+  is  the  temperature  of  the  deposit  just 
above  the  heated  rod  at  r  -  ru  and  T2~  is  the  temperature  of  the  deposit  just  below  the  fuel  at  r  =  r2. 
The  boundary  conditions  for  this  equation  are 


dj\ 

dr 


K, 


dIl 

dr 


at  r 


(ID 


and  the  condition  given  by  Equation  (7)  at  r  =  rv  In  Equation  (11),  the  subscript  3  denotes  the  fluid 
phase.  When  the  heat  flux  conditions  are  matched  at  r  =  rx  and  r  =  r2,  then 
TS=  Ti(r  =  ri)  =  T2(r  =  rx),  and  7*2~  =  T2{r  =  r2)  =  T3(r  =  r2). 

To  solve  Equation  (10),  the  thermal  conductivities  of  the  deposit  material,  Kd,  and  of  the  fuel,  Kf, 
are  required.  In  the  present  study,  we  use  JP-5  as  the  jet  fuel  and  the  properties  of  JP-5  as  a  function 
of  temperature  were  explicitly  employed  from  the  curve  fits  provided  to  us  by  the  sponsor.  Thus,  the 
thermal  conductivity  of  JP-5  in  (7 /mK  -  s)  is 

kj  =  A  +  BT  +  CT 2  (12) 

where  A  =  0.2503,  B  =  -5.6363  x  10-4,  and  C  =  4.7482 xlO-7.  We  further  assumed  that  the  thermal 
conductivity  of  the  deposit,  (/cd),  is  also  related  to  temperature  by  the  same  functional  relation  as  in 
Equation  (12),  but  since  the  deposit  thermal  conductivity  is  probably  much  lower  than  that  of  the 
fluid,  we  assume  Kd  =  Adnf  .where  Ad  <  1  is  a  constant  input  parameter.  Using  the  relation  between 
Kd  and  Kf  in  Equation  (10),  a  cubic  equation  is  obtained  which  can  be  written  in  the  nondimensional 
form  as 

^(7V>  -  7\3)  +  ^(r23  -  TS)  +  C*(T2  -  rx)  =  Cd  ln(I  +  h)  (13) 
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where  A*  =  AT? /Kdi,  B*  =  BTJk^,  and  C*  ~  C/Kd,  are  constants.  Here  h  =  (r2  -  rJ/rj  is  the  non- 
dimensional  thickness  of  the  deposit,  which  is  a  variable  in  this  problem.  Since  the  deposit  growth  is 
both  time-  and  space-dependent  in  the  present  formulation,  the  thickness  h  =  h(x,  t)  and  r2  =  r2(x,  t). 
In  Equation  (13)  we  have  neglected  the  superscripts  +  and  -  for  simplicity  and  the  nondimensional 
temperatures  Ti  and  T2  are  the  temperatures  at  7  =  1  and  r  =  1  +  h,  respectively.  At  this  stage,  there 
are  three  unknowns  in  Equation  (13),  i.e.,  7\,  T2,  and  Cd.  To  determine  these  unknowns,  we  combine 
Equation  (13)  with  the  boundary  conditions.  Equation  (7)  at  7  =  1  and  Equation  (11)  at  7  =  1  +  h. 
The  nondimensional  form  of  Equation  (7)  can  be  written  using  standard  one-sided  difference  as 

T1(7=l)  =  T1(F=l-Ar1)  +  alCd  (14) 

Similarly,  Equation  (11)  can  be  written  as 

T2(7  =  1  +h)  =  T2(7  =  1  +h  +  Ar2)  -  a2Cd  (15) 

where  the  differential  form  of  Equation  (10)  has  been  used  and 


[*/,  j  [(1 +£)*,] 


Here,  k,  and  Kf  are  the  nondimensional  thermal  conductivities  of  the  rod  and  the  fluid,  respectively. 
Also,  A rx  and  A r2  are  the  nondimensional  grid  spacing  in  the  rod  and  in  the  fluid,  respectively.  Note 
that  in  Equations  (14)  and  (15),  only  first-order  accurate  differencing  is  employed.  However,  a 
second-order  accurate  one-sided  difference  can  be  used  and  in  fact  was  employed  in  the  code.  Com¬ 
bining  Equations  (13)  through  (15),  we  obtain  a  cubic  equation  for  Cd.  The  real  root  of  this  equation 
is  the  solution  for  Cd.  Once  Cd  is  known.  Equations  (14)  and  (15)  can  be  used  to  determine  the  tem¬ 
peratures  at  r  =  1  and  7  =  1  +h.  Note  that  when  there  is  no  deposit  ( h  =  0),  then  Tx  =T2,  i.e,  the 
fuel  and  the  rod  are  at  the  same  temperature  at  their  interface.  However,  when  h  *  0,  the  cubic  equa¬ 
tion  must  be  solved  at  each  time  step  and  at  each  axial  location  to  match  the  gradients  and  thus  obtain 
the  temperature  of  the  rod  at  7  =  1  and  the  temperature  of  the  fuel  at  7  =  1  +h. 

An  advantage  of  this  formulation  is  that  it  allows  us  to  adjust  the  temperatures  of  the  rod  and  the 
fuel  depending  upon  the  heat  transfer  characteristics  of  the  deposit.  Thus,  when  Ad  <  1  is  given,  the 
deposit  will  inhibit  the  heat  transfer  from  the  rod  to  the  fuel  and  the  temperature  of  the  rod  at  7  -  1 
will  increase  while  the  temperature  of  the  fuel  at  7  -  1  +  A  will  decrease.  This  is  exactly  the  physical 
phenomena  observed  in  the  experiments  by  Marteney  (1988). 

3.3  Mass,  Momentum,  and  Heat  Transfer  in  the  Fluid  Phase 

The  fuel  deposit  will  form  due  to  the  degradation  reactions  that  occur  when  the  fuel  is  heated. 
Since  the  heat  source  is  a  metal  rod,  the  deposit  will  form  near  the  surface  of  the  rod  due  to  the  sur¬ 
face  catalytic  mechanism  This  catalytic  mechanism  is  very  complex  and  depends  on  various  factors, 
including  the  properties  of  the  fuel  and  the  surface  material  (e.g.,  Taylor,  1969b).  We  shall  discuss  this 
mechanism  in  a  later  section.  In  addition  to  this  mechanism  of  deposit  formation,  if  the  fuel  is 
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thermally  unstable  in  the  fluid  layer  near  the  high-temperature  surface,  then  the  fuel  can  undergo 
degradation  reactions,  resulting  in  deposit  formation  in  the  bulk  fuel  itself.  The  deposit  that  forms  in 
the  liquid  fuel  could  have  two  very  different  characteristics.  It  is  reasonable  to  assume  that  under 
some  conditions,  the  deposit  that  forms  due  to  chemical  reactions  in  the  bulk  fuel  will  be  liquid  pro¬ 
duct,  i.e.,  the  deposit  can  be  considered  a  liquid  precursor  which  may  eventually  undergo  a  phase 
change  and  condense  into  solid  form.  In  addition  to  this  liquid  deposit  precursor,  solid  particles  simi¬ 
lar  to  the  deposit  on  the  heated  wall  may  also  form  in  the  bulk  fuel.  Both  the  deposit  precursor  and 
the  solid  particles  can  migrate  to  the  wall  due  to  diffusion  processes  and  become  a  part  of  the  deposit 
on  the  surface.  A  reverse  migration  process  is  also  possible  due  to  the  fluid  motion  in  which  the  solid 
deposit  on  the  wall  is  entrained  into  the  flow  due  to  abrasion  and  flaking.  For  example,  in  wing  fuel 
tanks  containing  puddles  of  residual  liquid  hydrocarbon  and  hydrocarbon  vapor,  the  fuel  deposits 
have  a  tendency  to  flake  off,  contaminate  the  fuel,  and  cause  the  fuel  system  to  malfunction  (e.g., 
Taylor,  1969a). 

Thus,  to  model  the  fluid  phase  transport,  we  need  to  consider  the  general  mass,  momentum,  and 
heat  transfer  in  a  multicomponent  chemical  fluid  mixture  in  which  some  solid  particles  may  also  be 
present.  The  model  equations  must  be  a  generalized  version  of  a  two-phase  mixture  of  fluid  and  solid 
particles.  In  this  section,  we  will  develop  a  general  model  for  this  problem  which  allows  for  the  phase 
transition  of  the  liquid  deposit  precursor  into  solid  particles.  However,  since  this  is  a  complex  model, 
we  will  not  implement  it  in  this  phase.  After  deriving  the  general  model,  we  will  reduce  it  to  a  form 
which  is  simple  enough  to  be  considered  in  the  Phase  I  study  described  here. 

Before  deriving  the  general  governing  equations,  some  basic  assumptions  are  required  to  maintain 
the  tractability  of  the  equations.  Although  some  of  these  assumptions  are  quite  reasonable,  some  will 
have  to  be  reexamined  in  a  later  study  when  additional  experimental  data  is  available. 

The  primary  assumptions  used  to  derive  these  two-phase  flow  equations  are: 

1.  The  "bulk  fuel"  is  a  mixture  of  three  species:  the  fuel,  the  dissolved  oxygen,  and  the  liquid 
deposit  precursor. 

2.  Solid  deposits  may  be  present  in  the  "bulk  fuel"  due  to  phase  change  of  the  liquid  precursor, 
but  are  assumed  to  occupy  only  a  small  fraction  of  the  total  volume  of  the  bulk  fuel. 

3.  The  solid  deposit  in  the  "bulk  fuel"  is  mixed  well  with  the  fuel  so  that  the  solid  particles  can 
be  considered  a  pseudofluid,  thereby  allowing  the  derivation  of  the  solid  particle  transport 
equations  from  the  continuum  theory. 

4.  AH  solid  particles  in  the  "bulk  fuel”  are  of  the  same  size  and  can  be  approximated  by  a 
spherical  shape  of  radius  rp.  Thus,  the  drag  coefficient  of  all  the  particles  is  the  same  in  the 
fuel. 

Before  deriving  the  conservation  equations,  some  fundamental  definitions  must  be  given.  We  con¬ 
sider  a  two-phase  mixture  of  "bulk  fuel"  and  solid  particles.  It  must  be  noted,  however,  that  the  "bulk 
fuel"  itself  is  a  mixture  of  three  species  (fuel,  dissolved  oxygen,  and  liquid  deposit  precursor).  Dif¬ 
ferentiation  between  the  three  species  will  be  made  only  when  deriving  the  energy  transport  and  the 
bulk  fuel  species  transport  equations. 

We  define  a  solid  particle  in  the  bulk  fuel  with  a  mass  mp,  a  radius  rp,  and  a  specific  heat  C,,  all  of 
which  are  assumed  constant.  An  element  of  the  mixture  of  bulk  fuel  and  solid  particles  has  a  total 
mass  M  =  M}  +  Mp,  and  a  total  volume  V  =  V}  +  Vp,  where  the  subscripts  /  and  p  denote  the  bulk 
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fuel  and  solid  particles,  respectively.  If  the  total  number  of  solid  particles  in  the  bulk  fuel  is  given  by 
its  number  density  per  unit  volume,  tip,  then  the  total  volume  occupied  by  all  the  particles,  Vp,  is 

Vr  =  n  pV~*rp3  =  rtpVvp  (16) 

where  vp  is  the  volume  of  a  single  particle.  The  total  mass  occupied  by  all  the  solid  particles  is  given 
by 

Mp  =  mpnpV  (17) 

Thus,  the  species  density  of  the  solid  particles  is 


Pn>  = 


Hjl 

Vp 


08) 


For  a  given  fuel,  the  species  density  of  the  solid  particle  may  be  assumed  a  known  constant. 

For  a  two-phase  mixture,  we  need  to  know  the  partial  densities  of  the  pseudofluid  of  the  solid  par¬ 
ticles  and  the  fuel.  We  define  the  partial  density  of  the  particles  as 


-  Mp 

Pp  =  "TT  ~  rnpnp  =  Pn>vpnp 


(19) 


and  define  the  volume  fraction  of  the  solid  particles  in  the  two-phase  mixture.  A,  as 


A  = 


Yjl 

V 


npvp 


Pp_ 


(20) 


In  this  formulation,  A  is  a  principle  dependent  variable  representing  the  pseudofiuid  of  solid  particles 
in  the  bulk  fuel.  When  A  «  1  the  equations  can  be  further  reduced  to  the  form  representing  dilute 
mixtures  (Pai,  1977).  However,  in  the  present  formulation,  A  <  1  may  be  possible.  The  theory 
described  here  is  valid  for  finite  values  of  A. 

Similar  to  the  above  definitions,  we  define  the  species  density  of  the  bulk  fuel  as 


and  its  partial  density  as 


Pi  = 


MjL 

Vf 


Vf  V 


(1  -  A )p, 


(21) 


(22) 


With  these  basic  definitions  we  can  derive  the  conservation  equations  for  a  mixture  of  bulk  fuel  and 
pseudofiuid  of  solid  particles. 


TR-502/03-90 


10 


3.3.1  Equation  of  State 

Each  species  will  have  its  own  equation  of  state.  For  the  solid  particles  in  a  given  fuel,  we  may 
assume 

P, p  =  constant  (23) 

For  the  bulk  fuel,  the  equation  of  state  may  be  prescribed  as 


Pj  =  Pf(T j)  (24) 

where  Tf  is  the  temperature  of  the  bulk  fuel.  This  relationship.  Equation  (24),  may  be  assumed  given 
for  a  given  fuel.  For  example,  from  the  curve  fit  for  JP-5,  the  fuel  density  in  (kg  /m3)  is 

p ,  =  810.37  +  0.400477  -  1.3582  xl(T377a 


3.3.2  Conservation  of  Mass 

The  conservation  of  mass  for  both  phases  of  the  mixture  can  be  written  separately  as  follows. 


Solid  Phase 


dptpk  dPtpkUp. 

at  dxt  ~  °p 


(25) 


Here,  the  relation  =  pv A  has  been  used,  and  since  pv  is  assumed  a  constant  [Equation  (23)],  Equa¬ 
tion  (25)  is  an  equation  for  A.  Also,  uK  are  the  velocity  components  in  the  three  spatial  directions  of 
the  pseudofiuid  of  solid  particles  in  the  two-phase  mixture.  The  source  term  crp  indicates  that  there 
may  be  a  phase  transition  of  the  precursor  fuel-soluble  product  into  solid  particles.  If  op  =  >  0,  then 
the  precursor  fuel-soluble  deposit  will  convert  to  solid  particles.  At  present,  we  will  not  address  the 
functional  relationship  of  <rp. 


Bulk  Fuel  Phase 


dpf(  1-A)  dpj(  \-X)u,; 

dt  +  dx, 


=  -  cr„ 


(26) 


Here,  ufi  are  the  three  velocity  components  of  the  bulk  fuel.  The  mass  sink  in  the  bulk  fuel  continuity 
equation  implies  that  during  phase  transition  of  the  liquid  precursor  into  solid  deposit  there  will  be 
equivalent  loss  in  the  bulk  fuel  mass. 


Note  that  if  we  add  Equations  (25)  and  (26)  and  define  a  mixture  density  pM  as 

Pm  =  FP  +  Pi  =  kp*>  +  (1  -  k)p/ 
we  get  the  conservation  of  mass  of  the  mixture  as 

d  pM  dpMUMi 


dt 


dXi 


=  0 


where  the  mixture  velocity  u ^  is  given  as 


UU<=  \ppUK+  PfUf^/pu 


(27) 


(28) 


(29) 
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3.3.3  Conservation  of  Momentum 


The  conservation  of  momentum  for  both  phases  are  also  written  separately. 


Solid  Phase 


dXpvuK 

at 


qx  ^  Pp^ii  ryP]  —  F  Pi  F  p°i 


where  pp  =  X p  is  the  partial  pressure  of  the  solid  phase  and  p  is  the  pressure  in  the  mixture.  Also,  r,/ 
is  the  viscous  stress  tensor  of  the  pseudofluid,  which  by  definition  of  a  pseudofluid  may  be  modeled  in 
the  classical  form  as  a  function  of  the  viscosity  of  the  pseudofluid  (pp)  and  the  gradients  in  the 
pseudofluid  velocity  (wK).  The  forces  FK  and  F^.  are,  respectively,  the  i-th  component  of  the  solid- 
liquid  interaction  force  and  the  force  resulting  from  the  mass  source  ar 


Bulk  Fuel  Phase 


d{l-X)pJufl 


+  [o  +  pAi  -  V]  =  Fn  +  Fi°i 


Here,  p}  =  (1  -X)p  is  the  partial  pressure  of  the  bulk  fuel  and  is  the  viscous  stress  tensor  for  the 
bulk  fuel,  which  is  a  function  of  the  fuel  kinematic  viscosity  and  the  velocity  gradients.  F f.  is  the 
interaction  force  between  the  fluid  and  the  pseudofluid  of  particles,  and  F t„.  is  the  force  due  to  the 


mass  source  <r„ 


As  a  first  approximation,  the  viscous  stress  tensor  for  both  the  pseudofluid  of  solid  particles  and 
the  bulk  fuel  can  be  expressed  from  continuum  theory  as 


r  “  = 
7y 


where  n  stands  for  either  the  fluid  (/)  or  the  solid  phase  (p ),  *»«,  is  the  viscosity  of  the  nth  species,  and 
=  -2/3pni  is  the  bulk  viscosity. 


By  the  third  law,  we  have 


F,  =  -F„ 


F»  =  -F, 


The  interaction  force  (F {  or  Fp)  is  usually  a  very  complex  function  of  the  drag  coefficient  ( CD ),  the 
radius  of  the  particles,  the  number  density  (np),  and  the  difference  in  the  velocities  («p.  -  u}<).  The 
functional  form  of  this  force  will  have  to  be  prescribed  to  close  the  equations.  As  a  first  approxima¬ 
tion,  we  may  take  for  very  slow  flow  (i.e.,  for  low  Reynolds  number  flows),  the  Stokes  flow  approxi¬ 
mation  for  the  interaction  force  (Pai,  1977): 


FPi  =  6 *r„ii,np  [uf  ~  «;] 
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where  we  have  assumed  that  the  particle-particle  interaction  can  be  neglected  and  A  <  0. 1  applies,  so 
that  the  Stokes  form  for  drag  is  applicable.  When  0.1  <  A  <  1  occurs,  then  some  modification  to  the 
Stokes  law  will  have  to  be  made.  There  are  some  representative  expressions  for  the  interaction  force 
in  the  two-phase  flow  literature  which  could  be  used  here  (Pai,  1977). 

The  momentum  source,  f}  (or  FP(t\  is  a  function  of  the  mass  source  ap  when  there  is  phase 
change.  As  a  first  approximation,  we  may  consider  Fp<r  «  opi?p. 

3.3.4  Conservation  of  Energy 

Again,  the  conservation  of  energy  of  both  the  pseudofiuid  of  solid  particles  and  the  bulk  fuel  can 
be  written  separately. 


Solid  Phase 


We  first  define  the  total  energy  of  the  pseudofiuid  of  the  solid  particle  Ep  as  the  sum  of  its  internal 
energy  ep  and  its  kinetic  energy  as 

Ep  =  ep  +  jupj  (35) 

so  that  the  energy  conservation  equation  can  be  written  as 


dt  [^*p^p]  +  gx  [■^»pwpi-^p  ( uPjTiiP  dp,)  +  —  kt(F/  Tp)  +  £p  (36) 

where  ep  -  C,TP  is  the  internal  energy  of  the  pseudofiuid,  ep  is  the  energy  source  due  to  chemical 
reactions  or  heat  addition,  and  *Tis  an  interaction  coefficient  which  is  typically  a  function  of  rp,  np, 
and  kp.  Consistent  with  the  Stokes  law  approximation  used  in  the  interaction  force  for  the  momentum 
equation,  the  relation  kt  =  4*7ycp«p  can  be  used  (Pai,  1977).  Here,  np  is  the  thermal  conductivity  of 
the  solid  particles. 


In  Equation  (36),  qp.  is  the  heat  conduction  term  which  is  given  as 


dp,  = 


(37) 


Bulk  Fuel  Phase 

Similar  to  the  pseudofiuid  energy  equation,  an  energy  equation  for  the  bulk  fuel  can  be  written  as 
“  [(!  -*)/>,£/]  +  ~  [(1  -  A )P{u,£p  -  +  (1  -  A )u,;p\  =  kt{Tp-T,)  +  £/  (38) 

Note  that  the  interaction  term  appearing  on  the  right-hand-side  of  this  equation  is  the  negative  of  the 
interaction  term  in  the  pseudofiuid  energy  conservation  equation.  The  meaning  of  all  terms  in  this 
equation  are  similar  to  those  in  the  pseudofiuid  equation. 
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The  definition  of  the  total  energy  for  the  bulk  fuel,  Ef,  is  similar  to  that  for  the  pseudofluid  (i.e., 
Ef  =  ej  +  Uf?/ 2).  However,  since  the  bulk  fuel  is  considered  a  mixture  of  three  species,  some  further 
modification  to  the  internal  energy  of  the  bulk  fuel,  e} ,  must  be  made.  Thus,  when  there  are  three 
species  in  the  bulk  fuel,  the  basic  definition  of  internal  energy 

Plef  =  PfcvTf  (39) 

must  be  reinterpreted  as 

Ple!=  T,Pheh  (4°) 

*= 1 

where  k  =  1,3  represents  the  three  species  in  the  bulk  fuel  and  ejk  are  the  k-th  species  internal 
energy.  If  we  define  a  mass  fraction  Yk  =  p/k/pp  then  the  internal  energy  of  the  mixture  of  three 
species  becomes 

e/  =  i,YkCnT,  (41) 

*= i 

where  Cn  is  the  specific  heat  of  the  k-th  species. 

Equations  (23-26),  (30),  (31),  (36),  and  (38)  are  the  general  conservation  equations;  provided  that 
all  the  other  auxiliary  relations  are  defined,  these  equations  can  be  solved  for  the  dependent  variables 
A,  np.,  uf.,  pj,  Tp,  and  T f.  If  the  bulk  fuel  is  assumed  to  be  a  mixture  of  three  species,  then  additional 
equations  for  the  species  conservation  in  the  bulk  fuel,  Yk,  must  be  solved.  This  is  described  in  the 
next  section. 

Note  that  in  the  present  formulation,  we  have  assumed  that  the  bulk  fuel  is  laminar.  However,  in 
high  Reynolds  number  flows,  the  fluid  may  be  turbulent.  In  these  cases,  the  laminar  transport  proper¬ 
ties  (i.e.,  p.j,  kj,  etc.,  in  these  equations)  can  be  replaced  by  the  sum  of  the  laminar  and  turbulent  tran¬ 
sport  properties.  The  turbulent  properties  can  then  be  modeled  by  a  conventional  turbulence  model 
such  as  the  k  -e  model.  This  would  increase  the  number  of  equations  to  be  solved  in  the  bulk  fuel  and 
increase  the  computational  overhead.  For  fuel  flows  in  which  turbulent  effects  cannot  be  ignored, 
however,  a  turbulence  model  must  be  included.  This  could  be  considered  in  the  next  phase. 


3.3.5  Reduction  of  the  Two-Phase  Equations  for  the  Phase  I  Model 


The  equations  derived  above  are  very  complex  and  the  full  treatment  of  the  two-phase  problem  is 
beyond  the  scope  of  the  Phase  I  study.  Therefore,  we  consider  a  much  simpler  model  in  the  first 
phase.  The  first  basic  assumption  is  to  assume  that  there  are  no  solid  particles  in  the  bulk  fuel,  i.e., 
A  =  0.  This  immediately  reduces  the  system  of  equations  by  half.  The  next  major  simplification  is  to 
assume  that  the  bulk  fuel  is  being  transported  by  a  prescribed  axial  velocity,  «/,  so  that  for  an  axisym- 
metric  configuration,  uf.  =  («/,0).  Further,  we  assume  uf  =  uf(r).  Using  these  assumptions  and 
further  relating  the  bulk  fuel  density  to  the  bulk  fuel  temperature,  i.e.,  p{  =  Pf(Tf)  by  a  prescribed 
relation  (from  the  curve  fits  for  JP-5),  we  can  decouple  the  energy  transport  equation  from  the  mass 
and  momentum  equations.  Therefore,  the  only  equation  that  needs  to  be  solved  in  the  bulk  fuel  phase 
is  the  energy  equation  for  the  bulk  fuel.  Using  these  assumptions  and  assuming  there  are  no  sources 
in  the  fluid,  the  conservation  equation  for  the  bulk  fuel  energy  can  be  reduced  to 


dQ[, 

dXi 


(42) 
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If  we  further  define  the  total  enthalpy  as  Hf  =  Ef  +  p/p /,  Equation  (43)  can  be  rewritten  as 


i<""'  -  '> +  -k,  w**  -  -i <43> 

Since  we  are  prescribing  the  fluid  velocity  at  present  as  uSi  =  («(r),0),  and  assuming  that  the  pressure 
is  constant  in  the  JFTOT,  we  can  further  reduce  this  equation  for  axisymmetric  flows  to 


dh 

,  3hf } 

a 

dT,' 

_  i  a 

dT f 

,  l  a 

dut 

pf 

at 

+  u— L 
dx 

dx 

*'17 

+  7  3T 

r*'l 7 

r/i/U/17 

(44) 


where  P/(,Tf)  is  given  and  the  static  enthalpy  of  the  fuel  is  related  to  the  total  enthalpy  by  the  relation 
H,  =  hf  +  Uf2/2.  The  static  enthalpy  is  related  to  the  temperature  of  the  fuel  by  the  relation, 
hf  =  Cp/T j.  If  the  bulk  fuel  is  considered  to  be  a  mixture  of  three  species,  however,  this  definition  is 
modified  to 


h,  =  (45) 

*=i 

where  Cn  is  k-th  species  specific  heat  at  constant  pressure. 

Using  Equation  (45)  in  Equation  (43),  we  finally  obtain  an  equation  for  the  temperature  of  the  fuel 
T j.  The  only  unknown  in  this  equation  is  the  species  mass  fraction,  Yk ,  which  will  have  to  be  deter¬ 
mined.  This  is  described  in  more  detail  in  the  next  section. 

Equation  (43)  is  solved  subject  to  proper  boundary  conditions.  These  are: 


K, 


*13 

n 

p 

II 

10 

(46) 

o 

II 

{5 1* 

at  r  =  r3 

(47) 

II 

at  x  =  0 

00 

s-/ 

0 

II 

at  x  =  l 

(49) 

where  T j  and  Td  are  the  temperatures  of  the  bulk  fuel  and  of  the  solid  deposit,  respectively. 
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3.4  Species  Transport  in  the  Fluid  Phase 


3.4.1  Conservation  of  Species 

Since  the  bulk  fuel  is  assumed  to  be  a  mixture  of  three  species,  species  transport  for  each  species 
must  be  modeled.  We  consider  the  situation  were  the  formation  of  the  liquid  deposit  precursor  is  due 
to  the  reaction  between  the  fuel  and  the  dissolved  oxygen,  as  shown  by 


Fuel  +  Oxygen  — *  Deposit  Precursor 


(50) 


at  a  rate  that  is  currently  modeled  by  a  simple  one-step  Arrhenius  rate.  More  experimental  data  is 
needed  to  validate  this  assumption.  Therefore,  the  general  conservation  equation  for  the  k-th  species 
can  be  written  in  terms  of  the  concentration  of  the  k-th  species  as 

■^r +  ir*  h +M'JCJ ) =  Pc*  ~ Lc*  (51) 


where  [Ct]  is  the  concentration  of  the  k-th  species  in  the  bulk  fuel  and  utki  is  the  diffusion  flux  velo¬ 
city  of  the  k-th  species  in  the  i-th  direction.  The  diffusion  velocity  is  usually  given  by  Fick’s  law. 
The  last  two  terms  are,  respectively,  the  production  and  destruction  terms  for  the  k-th  species.  For 
the  present  discussion,  we  identify  k  =  1  with  the  fuel  (subscript  fu),  k  =  2  with  the  dissolved  oxygen 
(subscript  O'),  and  k  =  3  with  the  deposit  precursor  (subscript  pre).  The  use  of  a  simple  Arrhenius 
rate  for  the  production  and  destruction  terms  is  a  simplification  of  the  complex  chemical  kinetics 
occurring  in  the  conversion  of  the  fluid  to  the  deposit  precursor.  As  a  first  approximation,  the  pro¬ 
duction  of  the  deposit  precursor  in  the  bulk  fuel  can  be  modeled  as 


P 


A^txv 


[CfJa[C0]b 


(52) 


where  [Ct J  and  [CD]  are  the  concentration  of  the  fuel  and  molecular  oxygen,  respectively,  and  the 
exponents  a  and  b  are  at  present  taken  to  be  unity.  The  collision  frequency  A ^  and  the  activation 
energy  E ^  need  to  be  modeled  for  each  fuel.  This  model  is  similar  to  the  one  proposed  by  Krazinski 
et  al.  (1989).  The  species  molar  concentration,  [C/J,  can  be  rewritten  in  terms  of  the  mass  fraction  by 
using  the  relation  [C/J  =  PjiYfJ/(MjJ,  where  Mf%  is  the  molecular  weight  of  the  fuel  species.  Thus, 
Equation  (52)  can  be  rewritten  in  terms  of  the  k-th  species  mass  fraction  as 


Pj 


V(p,DkVYk)  +  Mk(PYt  -LYt) 


(53) 


Here,  Fick’s  law  tf  species  diffusion,  i.e.,  u/ti  =  -(Dk/Yk)VYk,  has  been  employed  and  the  mixture 
mass  conservation  equation  has  been  used.  Also,  Dk  is  the  diffusion  coefficient  for  the  k-th  species  in 
the  bulk  fuel.  To  solve  this  equation,  the  diffusion  coefficient  must  be  known  for  all  the  species  in  the 
bulk  fuel.  Currently  available  data  on  Dk  for  the  fuel  is  very  limited,  and  no  such  data  is  available  for 
the  deposit  precursor.  Specially  designed  experiments  will  be  required  to  determine  the  diffusion 
coefficients. 
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The  production  of  the  deposit  precursor  is  at  the  expense  of  the  fuel  and  oxygen  in  the  bulk  phase. 
Thus,  the  destruction  of  the  fuel  and  oxygen  species  in  the  bulk  fuel  is  directly  proportional  to  PY 
such  that 


and 


PY„  = 

~al>PYf„ 

(54) 

Py0  = 

~aoPY„. 

(55) 

where  <*/»and  acare  coefficients  that  are  functions  of  the  molecular  weights  Mk  and  the  exponents  in 
the  reaction  rate  (a  and  b)  in  Equation  (52). 

In  addition  to  the  production  of  the  deposit  precursor  in  the  bulk  fuel,  there  will  be  a  loss  of  the 
precursor  when  it  is  (a)  transported  to  the  wall  and  undergoes  phase  transition  into  solid  deposits,  (b) 
directly  converted  into  solid  deposits  in  the  bulk  fuel  itself  (pseudofluid  mass  source  <7P),  and  (c)  when 
the  liquid  deposit  precursor  is  converted  into  a  fuel-soluble  product  that  no  longer  has  the  tendency  to 
form  deposit  and  is  lost  from  the  deposition  process.  At  present,  due  to  lack  of  detailed  understand¬ 
ing  of  these  loss  mechanisms,  the  above  three  mechanisms  are  modeled  by  a  single  loss  term,  Lyrrc, 
which  represents  the  overall  loss  of  the  deposit  precursor  from  the  bulk  fuel.  Similar  to  the  earlier 
study  (Krazinski  et  al.,  1989),  we  describe  this  loss  term  by  another  Arrhenius  term  with  a  different 
collision  frequency  B ^  and  activation  energy  EL  such  that 


L 


Yprt 


(56) 


Again,  B ^  and  E L  both  have  to  be  calibrated  for  a  given  fuel,  and  c  is  the  exponent  which  for  simpli¬ 
city  we  choose  at  present  to  be  unity.  Note  that  the  loss  term  is  only  present  for  the  deposit  precursor, 
since  there  is  no  such  mechanism  for  the  fuel  and  the  dissolved  oxygen.  Thus, 


lYj.  ~  ly0  =  0 


(57) 


3.4.2  Reduction  of  the  Model  for  the  Phase  I  Study 

Since  we  are  at  present  studying  a  reduced  model,  we  again  make  some  simplifications.  We  first 
assume  that  the  fuel  and  dissolved  oxygen  can  be  lumped  together  into  a  single  species  (species  fO)  so 
that  only  two  species  have  to  be  modeled  rather  than  three.  By  mass  conservation,  y,o  +  ^  =  i; 
therefore,  only  one  species  transport  equation  has  to  be  solved.  We  therefore  solve  the  conservation  of 
the  deposit  precursor  species  (Y^  =  Y)  with  the  axisymmetric  assumption  so  that  the  conservation 
equation  for  Y  can  be  written  as 


dY 

dY 

d 

_  dY 

L  1  d 

dY 

Pi 

dt 

+  U'te 

dx 

"D’,ar 

+  7* 

<yDr— 

+  Py~  Ly 


(58) 
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where  the  production  and  loss  terms  are  given  by 


and 


Py  =  p,(\-Y)AYcxp 


Ly  =  PfYByeXp 


(59) 


(60) 


Here,  we  have  assumed  that  the  molecular  weights  of  the  fuel+dissolved  oxygen  and  the  deposit  pre¬ 
cursor  are  the  same.  As  noted  above,  the  solution  of  this  equation  requires  that  the  parameters  A  Y 
and  By  as  well  as  the  activation  energies  £Pand  EL  must  be  calibrated.  The  calibration  process  will 
again  depend  on  the  availability  of  detailed  experimental  data.  Equation  (58)  is  solved  subject  to 
proper  boundary  conditions. 


The  boundary  conditions  for  the  liquid  deposit  precursor  in  the  fluid  phase  are 


dY  n 

IT  =  0 

at 

r  =  r3 

o 

II 

at 

x  =  0 

x  =  / 

—  -  0 

at 

dx 

(61) 


The  boundary  condition  at  r  -  r2  is  very  important.  As  mentioned  earlier,  this  boundary  condition 
will  determine  the  diffusion/entrainment  process  occurring  at  the  edge  of  the  deposit  layer. 

The  current  lack  of  understanding  of  the  conversion  process  and  of  the  mechanism  that  affects  the 
entrainment  process  precludes  the  development  of  a  complex  boundary  condition  at  the  deposit  edge. 
Most  likely,  the  entrainment  of  the  liquid  deposit,  and  possibly  the  solid  deposit,  back  into  the  flowing 
fuel  would  be  in  some  manner  related  to  the  turbulence  levels  at  the  deposit  edge  and  the  surface  tur¬ 
bulent  shear  stresses.  Since  at  present  we  are  not  solving  the  turbulent  process,  we  will  not  attempt  to 
formulate  any  entrainment  process.  As  a  first  approximation,  we  may  assume  that  almost  all  the  pre¬ 
cursor  deposit  that  reaches  the  edge  of  the  solid  deposit  will  convert  in  some  manner  to  the  solid 
phase.  Thus,  we  must  determine  the  value  of  Y  at  r  =  r2.  Making  use  of  the  governing  form  of  the 
conservation  equation,  we  assume  that  near  the  edge  of  the  solid  deposit  layer,  the  convective  motion 
due  to  the  velocity  and  the  diffusion  process  can  be  considered  negligible  so  that  at  r  =  r2  we  may 
obtain  the  boundary  condition 

P, -f1  =  Py-Ly  (62) 

Thus,  we  solve  this  time-dependent  equation  at  the  solid  deposit  edge  with  the  initial  condition  at 
t  =  0  and  Y(r  =  r2)  =  0.  The  solution  of  this  equation  then  gives  the  required  boundary  condition  at 
the  edge  of  the  solid  deposit  on  the  wall. 
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3.5  Formation  of  Deposit  Due  to  Fuel  Degradation  Reactions 


The  fuel  deposit  will  form  due  to  the  degradation  reactions  that  occur  when  the  fuel  is  heated.  As 
noted  earlier,  there  are  three  possible  mechanisms  that  will  contribute  to  the  growth  of  the  deposit 
layer  on  the  heated  wall:  (a)  the  surface  catalytic  reaction,  (b)  the  transport  and  phase  change  of  the 
deposit  precursor  formed  in  the  bulk  fuel,  and  (c)  the  transport  and  impaction  of  the  solid  deposit 
formed  in  the  bulk  fuel.  The  net  result  of  these  mechanisms  will  be  the  growth  of  the  deposit  on  the 
metal  surface.  In  the  present  formulation,  we  define  the  thickness  of  the  solid  deposit  on  the  wall  as 
the  basic  parameter.  If  pd  is  the  density  of  the  deposit,  then  pdh  =  pj(r2  -  rj  will  be  the  net  deposit 
formed  per  unit  area  on  the  heated  metal  surface. 

The  final  result  of  the  fuel  degradation  process  is  the  formation  of  a  thin  layer  of  deposit  on  the 
metal  surface  of  the  heating  rod.  The  dimensional  thickness  of  the  deposit,  h  =  r7  -  rx  is  a  variable  in 
the  deposit  formation  process  and  can  be  modeled  by 


dh_ 

dt 


-  /4„exp 


dh  ' 

4_ 

dh 

dt 

I 

Y 

dt 

at  r  =  r2  (63) 


where  Ta  is  the  wall  temperature  and  the  first  term  models  the  surface  catalytic  conversion  by  an 
activation  energy  £„.  Note  that  Aw  has  the  unit  of  velocity.  Both  Av  and  £„  must  be  calibrated  for  a 
given  fuel.  The  modeling  of  the  surface  catalytic  process  by  the  single  Arrhenius  rate  equation  is  a 
major  simplification  of  the  complex  degradation  process.  At  present  it  is  not  clear  how  the  surface 
catalytic  process  will  degrade  the  fuel  and  result  in  the  formation  of  the  solid  deposit.  Due  to  lack  of 
proper  experimental  data  on  this  process,  the  current  model  uses  a  simplified  process  to  simulate  the 
effect  of  the  hot  wall  on  the  deposit  formation  rate.  Furthermore,  since  the  wall  temperature  is  nearly 
the  same  (i.e.,  a  weak  function  of  time),  the  calculation  of  the  deposit  growth  using  the  surface  cata¬ 
lytic  reaction  [assuming  the  other  two  terms  in  Equation  (63)  are  zero]  would  result  in  a  linear  growth 
of  the  deposit  with  time.  However,  experimental  data  (Marteney,  1988;  Marteney  and  Spadaccini, 
1986)  indicate  that  the  overall  deposition  rate  on  the  wall  is  nonlinear  with  time  and  in  fact  suggests 
that  the  growth  rate  increases  with  time.  This  nonlinearity  may  be  due  to  a  more  complex  catalytic 
reaction  mechanism  at  the  wall.  Again,  lack  of  experimental  data  makes  it  difficult  to  determine  if 
this  is  the  mechanism. 


The  second  term  models  the  overall  deposition  due  to  the  fluid  transport  of  the  deposit  precursor 
to  the  wall.  As  a  first  approximation,  we  may  assume  that  the  growth  of  the  deposit  near  the  wall  will 
be  in  some  manner  proportional  to  the  loss  of  the  deposit  precursor  from  the  bulk  fuel,  i.e.,  dhYoc  LY, 
where  L  Y  is  given  by  Equation  (60). 


Neglecting  for  the  moment  the  third  mechanism,  the  total  deposition  process  can  then  be  modeled 
by 


dh 

-  =  ^.exp 


RT. 


+  Bw  AZ.  Y 


(64) 


Here,  Bw  is  a  dimensionless  number  that  may  or  may  not  be  a  constant,  and  A  is  a  characteristic 
length  scale  within  which  the  phase  change  from  the  liquid  precursor  to  the  solid  deposit  would  occur. 
For  lack  of  data,  we  assumed  that  this  length  scale  is  equal  to  the  grid  scale  Ar2  near  the  edge  of  the 
deposit.  Computations  using  this  model  indicated  that  the  resulting  growth  is  not  much  different  from 
the  linear  growth  predicted  by  the  first  term  alone.  Using  the  JFTOT  data  (UTC,  1988),  a  series  of 
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computations  was  carried  out  to  see  how  the  nonlinearity  could  be  modeled  in  the  growth  mechanism. 
It  was  then  determined  that  the  nonlinear  growth  mechanism  can  be  fitted  to  the  experimental  data  by 
three  possible  modifications  to  the  above  equation.  The  first  modification  would  be  to  rewrite  Equa¬ 
tion  (64)  as 


dh 

dt 


(Au  +  am/)ex  p 


+  BwALy 


(65) 


where  is  another  adjustable  constant  with  the  units  of  acceleration.  However,  this  modified  model 
appeared  to  have  no  physical  significance  other  than  the  possibility  that  the  surface  catalytic  reaction 
changes  with  time.  Since  there  is  no  data  to  support  this  mechanism,  it  was  eventually  discarded. 

It  was  then  decided  to  look  for  the  nonlinearity  by  considering  the  modification  to  as  a  function  of 
the  deposit  thickness  itself.  This  led  to  the  second  modified  model: 


dh 

dt 


1  + 


exp  - 


RT„ 


+  B  w  AI .  y 


(66) 


Here,  aWl  is  again  a  constant  with  the  unit  of  velocity;  fitting  the  model  to  the  data  resulted  in  a  value 
of  aWl  «  2 Aw.  The  term  h*  is  a  reference  value  for  the  deposit  thickness  which  is  approximated  as 
follows:  As  the  deposit  grows,  there  may  be  a  critical  thickness  beyond  which  the  deposit  growth  will 
affect  the  heat  transfer  process.  At  this  limit,  we  may  approximate  the  heat  balance  at  r  =  r2  as 

ATd  AT, 

Kd  «  k,  — — 

n 


where  AT d  and  AT,  are,  respectively,  the  temperature  difference  in  the  solid  deposit  and  in  the  bulk 
fuel.  Here,  is  the  thickness  of  the  thermal  boundary  layer  above  the  deposit.  This  then  gives 

h  *  «  —  6th 

K, 


AT d 


AT. 


(67) 


If  we  further  assume  AT d  «  AT,,  Equation  (67)  then  gives  a  rough  estimate  for  h*.  Note  that  the 
assumption  Kd  =  Ad  k.,,  where  Ad  <  1,  implies  h  *  «  A Using  h  =  h  *  in  Equation  (66),  we  find  the 
characteristic  time  scale  to  reach  the  critical  thickness  as 


r 

' 

' 

A„[ 

I  +  cwij  exp 

RTW 

+  BkAL y 

The  model  equation  (66)  appears  to  have  some  significance  in  that  the  nonlinear  correction  to  the 
surface  catalytic  reaction  may  be  a  result  of  the  fact  that  the  surface  reaction  mechanism  may  be 
modified  by  the  growth  of  solid  deposit  on  the  wall.  Comparing  Equations  (65)  and  (66),  we  note  that 
both  these  models  can  be  considered  equivalent  if  we  assume  that  for  the  nonlinear  growth,  the  depo¬ 
sit  thickness  hat2.  However,  since  there  is  no  experimental  data  to  validate  this  possible  mechanism, 
we  were  unable  to  formulate  a  physical  process  that  would  be  consistent  with  the  Vh  variation  for  the 
nonlinear  process. 
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After  some  additional  consideration,  we  determined  that  the  nonlinear  term  could  also  be  modeled 
by  a  third  mechanism,  which  can  be  expressed  by 


dh  . 
~dt  = 


RTm 


+  Bm  A 


(68) 


This  mechanism  assumes  that  the  nonlinear  growth  of  the  deposit  is  in  some  way  related  to  the  deposi¬ 
tion  process  associated  with  the  phase  change  of  the  liquid  precursor  in  the  bulk  fuel  into  solid  depo¬ 
sit  near  the  wall.  After  some  theoretical  considerations,  we  eventually  formulated  a  possible  explana¬ 
tion  for  this  mechanism,  which  is  discussed  below. 


There  is  experimental  evidence  that  the  deposition  process  is  not  uniform  in  space  and  in  fact 
could  result  in  a  complex  shape  of  the  deposit  surface.  There  is  also  some  evidence  that  the  solid 
deposit  is  initially  very  porous.  Figure  2  shows  the  local  surface  of  a  deposit,  indicating  that  there 
may  be  regions  with  "holes"  into  which  the  bulk  fuel  and  the  liquid  deposit  precursor  may  be 
entrained.  In  the  regions  where  the  thickness  is  small,  the  bulk  fuel  will  be  in  contact  with  the  hot 
wall  and  the  surface  reaction  mechanism  will  result  in  a  locally  larger  deposition  rate.  This  mechan¬ 
ism  may  eventually  (after  a  long  time)  equilibrate  and  a  uniformly  thick  solid  deposit  would  form  on 
the  heated  surface.  However,  the  initial  surface  undulations  may  be  a  consequence  of  a  process  in 
which  the  liquid  deposit  precursor  in  the  bulk  fuel  undergoes  a  phase  change  into  solid  deposit  near 
the  wall  and  randomly  falls  and  adheres  to  the  metal  surface.  Thus,  in  general,  at  any  given  time  we 
have  an  "average"  deposit  thickness  defined  as  { h(x ,  t))  and  the  variance  of  h  given  by  »/h(x,  t). 

We  then  assume  that 


m. 

dt 


«  Vh 


(69) 


This  equation  implies  that  the  growth  of  the  average  deposit  thickness  with  time  will  increase  with  the 
increase  in  the  variance  or  the  extent  of  nonuniformity.  Since  we  are  assuming  that  the  nonuniformity 
of  the  deposit  is  a  result  of  some  random  process,  we  can  assume  that  this  random  process  can  be 
modeled  as  a  Gaussian  distribution.  If  P  is  the  probability  density  function  that  defines  the  probabil¬ 
ity  of  finding  the  thickness  within  the  range  h  and  h  +  dh,  then  by  definition 


(h)  =  jPhdh 

1*  =  Ip[h~  (h)}2dh 


(70) 


We  now  need 
assume  that 


to  define  the  probability  density  function  in  terms  of  the  Gaussian  distribution.  We 


P 


(71) 


where  G  is  the  Gaussian  probability  distribution.  By  definition,  the  Gaussian  distribution  satisfies 


and 


jGdh  =  1 


(72) 


(h)  =  jGhdh 

Ik7  =  !G\h-{h)}*dh 


(73) 
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It  can  be  shown  that  this  definition  satisfies  Equation  (70)  by  considering  the  time  derivative  of  Equa¬ 
tion  (73)  as 

^  -  /f*»  <74> 

Using  the  definition  of  P  from  Equation  (71),  we  can  easily  show  that  Equation  (74)  is  an  identity,  i.e., 
the  right-hand  side  becomes  the  same  as  the  left-hand  side.  If  we  follow  the  same  procedure  and  com¬ 
pute  the  first  time-derivative  of  the  variance,  »jka,  we  obtain 

<75> 

Integrating  Equation  (75)  results  in 

nS  =  ?{h) 

or 

fl  h  =  (76) 

Thus,  by  Equation  (68),  we  obtain 

*  Vh  a  P(h)2  (77) 

Here,  p  is  a  constant  of  the  order  of  unity.  Thus,  the  nonlinear  term  in  Equation  (67)  can  be  related 
to  the  random  deposition  process  near  the  wall  by  virtue  of  Equation  (77). 

Finally,  in  Equation  (63),  we  had  identified  a  third  possible  mechanism  for  the  growth  of  the  solid 
deposit  layer  on  the  wall.  This  mechanism  is  related  to  the  solid  particles  formed  in  the  bulk  fuel 
when  the  liquid  precursor  condenses  into  solid  particles.  These  solid  particles  will  then  migrate 
towards  the  wall  and,  upon  reaching  the  wall,  will  begin  to  adhere  to  it.  Note  that  this  mechanism  is 
different  than  the  mechanism  for  ( dh/dt )  discussed  earlier.  In  that  model.  Equation  (68),  the  addi¬ 
tional  solid  deposit  formed  on  the  wall  is  the  result  of  the  phase  change  of  the  liquid  precursor  into 
solid  deposit  near  the  deposit  surface.  Diffusion  through  the  bulk  fuel  is  not  important  for  this 
mechanism.  For  the  third  mechanism,  however,  the  migration  of  the  solid  deposit  formed  in  the  bulk 
fuel  towards  the  wall  is  important.  This  migration  and  subsequent  adhesion  process  needs  to  be 
modeled  because  it  will  contribute  to  the  overall  deposition  rate.  Here  we  use  as  a  first  approximation 
the  approach  used  in  aerosol  deposition  theories  (Davies,  1987)  and  formulate  the  third  mechanism  as 
a  function  of  a  diffusion  coefficient  and  the  normal  gradient  of  the  volume  fraction  of  the  solid  parti¬ 
cles  in  the  bulk  fuel.  Thus, 

=  C^-g-  (78) 

.  D 

where  is  the  diffusion  coefficient  for  the  solid  particles,  Ci  is  an  adjustable  constant  of  the  order 
of  unity,  and  the  equation  is  valid  at  r  =  r2.  From  aerosol  deposition  theory,  we  estimate  the  diffusion 
coefficient  in  terms  of  the  vertical  velocity  fluctuation  and  a  characteristic  length  scale.  Thus, 

Dv  =  v;  A  (79) 
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(80) 


where  the  vertical  velocity  fluctuation  near  the  deposit  edge  is  given  by 

v/  =  fiU*Sc * 

In  Equation  (80),  fx  is  a  constant  and  u  *  is  the  friction  velocity,  usually  given  as  u  *  =  Ip).  Typi¬ 
cally,  u*  is  a  small  fraction  of  the  free  stream  velocity,  ua,  Le.,  u*  «  fiu„  where  /2  «  0.01  -0.05. 
Here,  r„  =  pj(du/dr)  is  the  shear  stress  at  the  edge  of  the  deposit  (r  =  ra).  Also,  Sc  =  v/D^  is  the 
Schmidt  number  and  n  is  an  exponent.  The  value  of  fx  and  n  depends  in  part  on  the  state  of  the  wall, 
i.e.,  for  smooth  walls,  aerosol  deposition  data  indicates  /j  =  0.075  and  n  =  -2/3;  for  rough  walls, 
fx  -  0.08  and  rt  —  -1/2.  For  the  present  study,  the  rough  wall  assumption  is  probably  more  appropri¬ 
ate  and  thus  is  employed  here.  The  value  of  the  Schmidt  number  is  typically  in  the  range  of  106  for 
very  small  particles.  Therefore,  for  the  present  study  we  used  n  =  -2/3,  fx  =  0.08,  and  Sc  =  106. 

In  the  Phase  I  model,  we  had  assumed  that  there  is  no  solid  deposit  present  in  the  bulk  fuel.  Thus, 
Equation  (78)  cannot  be  used  for  the  third  mechanism  directly.  An  indirect  way  to  implement  this 
mechanism  is  to  rewrite  Equation  (78)  as 

%  =CwDy^r  at  r  =  r2  (81) 

D 

where  Cw  is  another  adjustable  constant  of  the  order  of  unity  and  DY  is  the  diffusion  coefficient  for 
bulk  fuel  precursor  near  the  deposit  edge,  which  is  assumed  to  be  the  same  as  D n,.  This  mechanism 
essentially  indicates  that  the  normal  diffusion  of  the  liquid  precursor  formed  in  the  bulk  fuel  towards 
the  wall  will  result  in  an  additional  growth  in  the  solid  deposit.  However,  the  liquid  precursor  can 
only  become  solid  deposit  by  undergoing  a  phase  change.  Thus,  if  the  temperature  near  the  solid 
deposit  surface  is  high,  the  liquid  deposit  precursor  will  not  solidify  and  instead  will  form  a  liquid 
boundary  layer  flowing  in  the  downstream  direction  on  top  of  the  solid  deposit.  This  liquid  layer  will 
Starr  to  solidify  only  if  the  temperature  drops  below  a  critical  point,  characteristically  the  condensa¬ 
tion  temperature.  Therefore,  we  make  an  additional  restriction  that  the  third  mechanism  contributing 
to  the  deposit  growth.  Equation  (81),  will  be  negligible  if  AT  =  Tw(x)  -  Tw(x  -  Ax)  >  0,  i.e.,  the 
temperature  is  increasing  in  the  downstream  direction,  and  ( dh/dt)D  *  0  only  if  AT  <  0,  i.e.,  the  fuel 
temperature  at  r  =  r2  starts  to  decrease  in  the  outflow  direction.  At  present,  this  modified  formula¬ 
tion,  although  reasonable,  cannot  be  fully  validated.  Therefore,  in  the  discussion  of  the  results  later  in 
this  report,  no  major  emphasis  is  placed  on  this  model. 

This  completes  the  formulation  of  the  system  of  differential  equations  that  governs  the  five  unit 
processes  involved  in  the  fuel  degradation  and  deposit  formation  process.  These  unit  processes  are 
coupled  through  the  boundary  conditions  at  the  interfaces.  If  detailed  data  were  available,  it  would 
be  worthwhile  to  evaluate  the  validity  of  the  full  model  formulated  here.  In  the  reduced  model  for  the 
Phase  I  investigation,  no  solid  particles  are  assumed  to  be  present  in  the  bulk  fuel.  The  velocity  field 
is  assumed  known,  and  the  transport  of  heat  is  a  passive  scalar  transport  process.  The  solid  phase  can 
be  added  in  a  subsequent  study  if  further  data  is  available  for  validation.  The  parameters  in  the 
present  model  are  the  thermal  conductivities  of  the  heating  element  and  the  deposited  material,  the 
activation  energy  and  rate  constants  for  the  surface  catalytic  reaction  and  the  fluid  phase  reaction, 
and  the  viscosity  and  thermal  conductivity  of  the  bulk  fuel.  The  operating  parameters  are  the  dimen¬ 
sions  of  the  apparatus,  the  flow  rate  of  the  fuel,  and  the  heat  flux  into  the  heating  element. 
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4.  THE  NUMERICAL  CODE 

The  reduced  mathematical  model  of  the  fuel  degradation  problem,  as  described  in  the  previous 
section,  was  implemented  in  a  numerical  code  to  carry  out  a  series  of  simulations.  For  the  eventual 
development  of  a  general  purpose  thermal  prediction  code,  we  undertook  the  task  of  developing  a 
completely  new  CFD  code  for  this  problem.  The  numerical  method  used  is  described  below. 

In  the  Phase  I  research,  the  fluid  was  assumed  to  be  JP-5  with  prescribed  transport  properties,  i.e., 
Pf(T),  cPj(T),  itf(T ),  n;(T),  and  with  a  known  flow  velocity.  Thus,  the  fluid  dynamics  is  entirely 
decoupled  from  the  heat  exchange  problem.  The  velocity  profile  for  a  flow  in  an  annular  tube  can  be 
obtained  exactly  from  classical  considerations  (Landau  and  Lifshitz,  1987)  as 


uf(r) 


_A£_ 
4  P,l 


(82) 


where  A p  is  the  imposed  pressure  gradient,  r3  is  the  outer  radius,  and  r2  is  the  inner  radius,  which  in 
the  present  case  of  deposit  growth  would  be  r2  =  r2(x,/).  Also,  l  is  the  length  of  the  tube  ana  p}  is  the 
kinematic  viscosity  of  the  fluid.  Strictly  speaking,  if  r2  =  r2(x,t)  for  the  cases  with  deposit  growth, 
then  Uj  will  vary  with  (*,/).  In  that  case,  the  velocity  profile  at  each  axial  location  and  at  each  time 
step  can  be  modified  to  take  this  into  account. 

The  heat  transfer  in  the  rod.  Equation  (2),  the  heat  transfer  across  the  deposit.  Equation  (9),  and 
the  heat  balance  in  the  fuel.  Equation  (44),  were  solved  using  both  an  explicit  time-marching  scheme 
and  an  implicit  ADI  scheme.  To  achieve  high  accuracy,  a  high-order  spatial  discretization  scheme 
was  developed.  Steady-state  temperature  fields  without  deposit  growth  were  calculated  using  both 
schemes. 


4.1  The  Explicit  Scheme 

The  explicit  scheme  consists  of  a  high-order  compact  scheme  (the  Pa'de  approximation  scheme) 
from  which  a  formally  fourth-order-accurate  scheme  can  be  obtained  by  using  a  compact  stencil.  If 
we  denote  60Uj  =  (ut+1  -  u,_i)/2Ax  as  the  second-order-accurate  central  difference  for  the  first  spa¬ 
tial  derivative  and  I?  =  (n,+1  -  2w,  +  «,•_ i)/Ax3  as  the  second-order-accurate  second  derivative,  then 
by  Taylor  series  we  have 

+  —t-  ~7  +  0( Ax*)  (83) 

1  dx  6  dx3 

Thus,  the  first  derivative  can  be  written  from  Equation  (83)  as 


du 

dx 


1+-2T*3 

o 


-l 


«.-  +  !  ~  Ui-l 

2  Ax 


(84) 


where  P  is  the  central  second  difference  operator.  This  discretization  is  on  a  uniform  grid.  Extension 
to  a  nonuniform  grid  is  straightforward.  To  solve  this  equation,  a  tridiagonal  system  must  be  inverted. 
Also,  since  specific  boundary  conditions  are  prescribed  in  the  present  case,  special  consideration  must 
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be  given  to  the  discretization  at  the  boundaries.  Instead  of  the  central  differences,  we  employ  one¬ 
sided  second-order  difference  near  the  boundaries  so  that 


-3 u0  +  4ui  -  u2 
2Ax 


du0 

dx 


^  +0(Ax4) 


(85) 


This  spatial  discretization  was  combined  with  the  Runge-Kutta  time-stepping  scheme  to  produce  a 
highly  accurate  unsteady  scheme.  The  scheme  developed  in  this  study  was  written  so  that  spatial  and 
temporal  accuracy  of  up  to  fourth-order  could  be  obtained.  A  second-order-accurate  scheme  was  also 
developed  in  this  study. 


Fourth-order  compact  schemes  can  have  problems  due  to  numerical  instability  in  flows  with  sharp 
gradients.  Both  second-  and  fourth-order  schemes  were  implemented  and  were  then  compared  using 
simple  test  problems.  For  one-dimensional  wave  motion  with  periodic  boundary  conditions,  it  was 
determined  that,  with  the  compact  scheme,  only  six  grid  points  are  required  to  resolve  the  wave  rea¬ 
sonably  without  significant  phase  errors,  whereas  the  second-order  scheme  required  a  much  higher 
resolution.  For  the  present  configuration  with  a  constant  heat  flux  into  an  initially  cold  rod  and  fuel, 
the  initial  sharp  temperature  gradient  (at  l  *  0)  caused  numerical  instability  when  the  fourth-order 
scheme  was  used.  When  the  initial  temperature  distribution  in  the  rod  and  in  the  bulk  fuel  was  more 
uniform,  the  compact  scheme  encountered  no  major  instability.  Thus,  it  appeared  that  the  best 
approach  would  be  to  begin  the  computation  using  the  second-order  scheme  and  then  switch  to  the 
higher-order  scheme  once  the  temperature  field  in  the  rod  and  in  the  fuel  has  developed.  For  the 
reduced  model  used  in  Phase  I,  both  the  second-  and  fourth-order  schemes  resulted  in  identical  steady 
states.  Since  the  computational  overhead  for  the  second-order  scheme  was  significantly  lower,  most  of 
the  steady-state  computations  employed  the  second-order  scheme. 

4.2  The  Implicit  Scheme 

The  basic  problem  with  the  explicit  scheme  is  that  the  stability  considerations  severely  limit  the 
time  step  for  the  heat  conduction  equation  to  very  small  steps,  which  leads  to  excessive  computational 
cost.  To  avoid  this  restriction  on  the  time  step,  we  implemented  an  implicit  scheme  based  on  the 
alternate  direction  implicit  (ADI)  scheme.  Such  schemes  are  well  known  and  have  been  extensively 
used  in  the  study  of  such  problems  where  the  time  step  restriction  can  be  quite  stringent.  We  will  not 
go  into  the  details  of  this  scheme  since  it  is  described  thoroughly  in  many  textbooks  (Anderson  et  al., 
1984).  Here  we  briefly  summarize  the  main  steps.  Consider  a  general  equation  of  the  form 


d<t> 

at 


+  b<t>x  +  c<f>n  +  d<f>r 


(86) 


where  a,  b ,  c,  and  d  need  not  be  constant  and  in  the  present  case  are  assumed  to  be  functions  of  4>. 
Here,  <f>  may  be  considered  to  represent  the  temperature,  T.  Also,  the  subscripts  indicate  differentia¬ 
tion  with  respect  to  x  or  r;  for  example,  the  subscript  xx  indicates  the  second  derivative  of  <i>  with 
respect  to  x.  Further,  we  will  denote  Dt  and  D,  as  the  finite  difference  approximation  to  the  x- 
derivatives  and  the  r-derivatives,  respectively.  Then,  using  the  Crank-Nicholson  scheme,  we  can 
write 


'  -  fD*  -  fD- 


I  +  —~(DX  +  Dr) 


Tn  +  0(A/2) 


(87) 
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Using  approximate  factorization,  we  obtain 


T* 


(88) 


In  the  present  case,  since  the  coefficients  in  Equation  (86)  are  functions  of  <f>,  we  make  the  simpli¬ 
fying  assumption  that  the  coefficients  and  the  boundaries  are  constant  for  one  time  step.  This  reduces 
the  time-accuracy  to  0(  A/)  but  will  not  affect  the  solution  if  only  steady-state  solutions  are  required. 


We  solve  this  equation  in  the  rod  using  the  Peaceman-Rachford  ADI  scheme  (Peaceman  and  Rach- 
ford,  1955)  as  follows 


7”* 


(89) 


•J*  »  +  l 


(90) 


A  similar  scheme  was  used  to  solve  for  the  temperature  in  the  fluid;  in  this  case  the  convection  term 
was  added  into  the  x-operator.  A  modification  of  this  type  was  used  by  Polezhaev  (1967).  Since  in 
each  time  step  the  coefficients  are  changing,  we  linearized  locally  and  assumed  that  the  coefficients 
and  boundary  conditions  remained  constant  during  each  time  step.  After  a  complete  time  step,  the 
coefficients  and  boundary  conditions  were  updated.  This  reduced  the  total  accuracy  to  0(A/,  Ax2). 

Since  there  is  a  convective  term  in  the  present  case  due  to  the  velocity  field,  this  scheme  will  not 
be  unconditionally  stable;  however,  the  time  step  restriction  with  dependence  on  Ax  for  the  implicit 
scheme  is  still  much  larger  than  the  time  step  possible  with  the  explicit  scheme  (Anderson  et  al., 
1984).  All  the  steady-state  results  presented  in  this  report  employed  the  implicit  scheme. 
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5.  RESULTS  AND  DISCUSSION 

In  this  section,  we  describe  the  results  of  the  study.  Before  describing  the  results,  it  is  important  to 
emphasize  that  there  were  quite  a  few  variables  that  had  to  be  calibrated  for  the  deposition  process. 
We  essentially  assumed  that  the  data  on  the  JFTOT  could  be  used  for  this  purpose.  However,  the  data 
itself  was  quite  sparse  and  there  were  many  parameters  that  were  still  unknown.  Therefore,  although 
the  model  developed  in  the  present  study  appears  to  have  the  proper  physics  to  simulate  the  deposi¬ 
tion  formation  and  the  growth  of  the  solid  deposit  on  the  wall,  the  full  evaluation  of  the  model  will 
require  additional  data  on  the  physical  processes  involved. 

The  JFTOT  data  used  here  (UTC,  1988)  for  calibration  had  three  distinct  features.  The  first  was 
that  the  measured  wall  temperature  on  the  JFTOT  was  not  uniform  and  in  fact  had  a  peak  at  a  fixed 
axial  location.  After  some  initial  computations,  it  became  quite  clear  that,  with  a  constant  heat  flux 
boundary  condition  at  the  centerline,  the  computed  temperature  of  the  wall  did  not  have  the  variation 
observed  in  the  JFTOT  data.  It  appears  that  the  peaked  temperature  distribution  on  the  wall  is  due  to 
the  three-dimensionality  of  the  flow  field  inside  the  JFTOT  (CRC,  1978;  Oh  et  al.,  1989)  and  to  the 
fact  that  the  fuel  outflow  region  is  cooled,  so  that  significant  heat  loss  occurs  at  the  outflow.  Since  the 
present  configuration  is  axisymmetric,  it  was  not  possible  to  predict  the  wall  temperature  distribution 
that  was  observed  in  the  JFTOT  experiments  by  using  the  constant  uniform  heat  flux  conditions.  It 
was  possible,  however,  to  mimic  the  3D  effects  by  modifying  the  heat  flux  condition.  This  modifica¬ 
tion  is  described  below.  For  the  calibration  studies,  however,  we  decided  to  use  the  curve-fitted  tem¬ 
perature  distribution  provided  by  the  sponsor  for  the  wall  temperature. 

The  second  feature  observed  in  the  JFTOT  data  (UTC,  1988)  was  that  as  the  duration  of  the  test 
increased,  the  deposition  rate  also  increased.  This  nonlinear  behavior  appears  to  be  in  agreement  with 
the  observations  of  Marteney  (1988).  As  mentioned  earlier,  this  nonlinear  growth  data  was  the  pri¬ 
mary  information  used  to  calibrate  the  model. 

Finally,  the  third  distinct  feature  observed  in  the  JFTOT  experiments  was  that  solid  deposit  distri¬ 
bution  was  biased  in  the  downstream  direction  such  that  more  solid  deposit  was  present  downstream 
of  the  peak  temperature  location.  Some  attempt  to  model  this  unusual  feature  was  attempted  and  is 
also  described  here. 

The  process  of  model  calibration  involved  first  obtaining  the  steady-state  temperature  distribution 
on  the  rod  and  in  the  fluid  using  the  assumption  that  no  deposit  was  forming.  This  essentially 
required  the  solution  of  Equation  (1)  and  (44)  subject  to  proper  boundary  conditions.  Since  the 
JFTOT  wall  temperature  distribution  could  not  be  reproduced  by  the  present  axisymmetric  model,  we 
obtained  the  steady-state  temperature  field  in  the  bulk  fuel  by  specifying  the  temperature  profile  near 
the  centerline  of  the  rod  as  the  boundary  condition  so  that  at  steady  state  the  wall  temperature  at 
r  =  rt  matches  the  experimental  curve  fit.  Once  the  steady-state  temperature  field  in  the  bulk  fuel 
was  determined,  we  solved  the  deposit  growth  equation  (68)  and  the  liquid  precursor  growth  equation 
(62)  with  the  temperature  field  given  by  the  steady-state  field.  The  effect  of  deposit  growth  on  the 
heat  transfer  process  was  taken  into  account  by  the  method  described  in  Section  3.2.  The  calibration 
of  the  data  for  the  given  temperature  field  was  then  carried  out.  The  results  of  these  studies  are  sum¬ 
marized  here. 

The  basic  JFTOT  configuration  used  in  this  study  had  the  following  dimensions  and  reference 
conditions: 

r1  -■  1.5625  mm,  r3  =  2.3125  mm,  /  =  60  cm 
ua  =  0.0107  m/s,  T„  =  293AT,  Re  «  10 
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All  other  reference  properties,  such  as  /j/co,  /c/oo,  and  pfoo,  are  computed  using  the  curve  fit  at  the 
inlet  temperature,  T0.  Notice  that  for  the  flow  speed  used  here,  the  Reynolds  number  is  very  low; 
therefore,  slow  flow  (Stokes)  approximation  could  be  used  to  neglect  the  convective  effects.  However, 
we  retain  the  convective  terms  in  our  formulation.  The  low  Reynolds  number  of  the  flow  implies  that 
laminar  flow  approximation  is  applicable  here. 

Steady-state  solutions  for  various  boundary  conditions  were  obtained  to  determine  the  capabilities 
of  the  explicit  and  implicit  numerical  schemes.  In  general,  the  solutions  obtained  with  the  two 
schemes  were  nearly  identical.  Since  the  implicit  scheme  is  more  appropriate  for  the  steady-state  cal¬ 
culation,  all  the  results  presented  here  employed  the  implicit  scheme.  Furthermore,  only  characteristic 
solutions  are  presented  here. 

Figure  3  shows  the  temperature  distribution  in  the  rod  and  in  the  fuel  for  the  case  in  which  a  con¬ 
stant  heat  flux  condition  was  used  near  the  centerline.  For  this  calculation,  64  grid  points  in  the  x- 
direction  and  42  grid  points  in  the  radial  direction  (10  points  in  the  rod,  32  points  in  the  fuel)  were 
used,  and  it  was  assumed  that  the  rod  is  insulated  at  x  =  0  and  x  =  /.  A  constant  heat  flux  of  400 
kW/m2  was  used  for  this  calculation.  Figure  3a  shows  the  temperature  surface  plotted  such  that  the 
fuel  inflow  temperature  was  normalized  to  zero.  The  constant  temperature  surface  shown  in  Figure  3a 
is  plotted  in  the  computational  plane  and  thus  is  not  to  scale.  Figure  3b  shows  the  temperature  con¬ 
tours  in  the  rod  and  in  the  fuel -again  plotted  in  the  computational  plane  for  clarity  (since  the  rod  is 
very  long  and  very  thin).  In  comparing  Figures  3a  and  3b,  we  note  that  there  is  a  temperature  boun¬ 
dary  layer  in  the  fluid  above  the  rod.  Since  the  fuel  at  the  inlet  is  cold  and  is  flowing  in  the  positive  x 
direction,  it  convects  the  heat  downstream  so  that  the  fuel  temperature  is  lower  near  the  inlet  and  the 
temperature  continuously  rises  towards  the  outflow.  In  Figure  3c,  the  temperature  profiles  are  plotted 
at  various  axial  locations  and  the  thermal  boundary  layer  is  more  apparent.  Peak  temperature  is 
observed  near  the  outflow. 

We  repeated  this  calculation  (not  shown  here)  with  the  assumption  that  there  are  heat  losses  at  the 
two  ends  of  the  rod.  However,  there  was  no  major  change  in  the  temperature  distribution  except  for  a 
decrease  in  temperature  near  the  outflow.  These  two  calculations  clearly  indicated  that  the  constant 
heat  flux  boundary  condition  will  not  result  in  a  temperature  distribution  as  observed  in  the  JFTOT 
data.  However,  the  temperature  field  distribution  obtained  in  these  calculations  were  similar  to  those 
observed  in  the  heated  tube  experiments  by  Marteney  (1988),  which  showed  that  the  temperature  con¬ 
tinues  to  rise  toward  the  outflow  region. 

These  calculations  were  then  repeated  with  the  heat  flux  modified  at  the  centerline  such  that  it 
had  a  parabolic  shape  with  a  maximum  heat  flux  at  the  midpoint  of  the  rod,  i.e.,  x  =  30  cm.  The  peak 
heat  flux  was  adjusted  so  that  the  peak  wall  temperature  was  near  the  observed  value  for  a  specified 
JFTOT  case  (UTC,  1988)  with  a  peak  temperature  of  287'C  (550'F).  Figures  4a  through  4c  show  the 
temperature  surface,  the  temperature  contour,  and  the  temperature  profiles,  respectively,  for  this  cal¬ 
culation.  Figure  4a  clearly  shows  that  there  is  a  peak  in  the  temperature  field  near  the  wall.  The  tem¬ 
perature  contours  shown  in  Figure  4b  indicate  that,  on  the  wall,  the  temperature  peak  is  not  at 
x  =  30  cm  where  the  input  was  a  maximum  but  rather  a  small  distance  downstream  of  the  midpoint 
location.  This  is  primarily  due  to  the  convective  effects  caused  by  the  fluid  flow  on  the  top  of  the  rod. 
This  downstream  shift  of  the  peak  location  is  thus  consistent  with  observations  in  the  JFTOT  experi¬ 
ments  (UTC,  1988).  The  thermal  boundary  layer  profiles  at  various  axial  locations  show  that  down¬ 
stream  of  the  peak  location,  the  thermal  boundary  layer  thickness  decreases  since  the  peak  values  at 


TR-502/03-90 


29 


Figure  3. 


a.  Steady-state  temperature  surface 


Steady-State  Temperature  Field  in  the  Heating  Rod  and  in  the  Fuel  with  a 
Constant  Heat  Flux  Boundary  and  Insulated  Ends 


TR-502/03-90 


30 


tjjix  Radius 

os*  sz-  oo  t  szm  os-i 


0 


X 


Figure  3. 

R-502/0V 


b.  Temperature  contours  in  the  rod  and  fuel.  Contour  intenal  is  20°C. 
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c.  Temperature  profiles  in  the  rod  and  fuel  at  various  axial  locations 


Steady-State  Temperature  Field  in  the  Heating  Rod  and  in  the  Fuel 
Constant  Heat  Flux  Boundary  and  Insulated  Ends  (Cont.) 
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a.  Steady-state  temperature  surface 


Figure  4.  Steady-State  Temperature  Field  in  the  Heating  Rod  and  in  the  Fuel  with  a 
Parabolic  Heat  Flux  Boundary  and  Insulated  Ends 
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the  wall  are  dropping.  With  the  modified  heat  flux  boundary  condition,  this  calculation  essentially 
mimicked  the  temperature  distribution  observed  in  the  JFTOT  data.  Since  the  actual  JFTOT  experi¬ 
mental  setup  was  cooled  at  the  downstream  end,  we  also  repeated  the  above-mentioned  calculations 
with  a  heat  loss  at  the  downstream  end  (not  shown  here).  No  major  change  in  the  solution  was 
observed  except  for  a  decrease  in  the  temperature  as  the  outflow  is  reached.  These  types  of  boundary 
heat  flux  modifications  are  somewhat  arbitrary  since  the  JFTOT  flow  field  is  three-dimensional  and 
the  current  study  is  axisymmetric. 

Subsequently,  we  decided  not  to  attempt  to  reproduce  the  JFTOT  wall  temperature  data  using  the 
constant  or  modified  heat  flux  boundary  conditions.  Instead,  we  used  the  curve  fit  provided  by  the 
sponsor  for  the  wall  temperature  distribution,  of  the  form 

Tj(x)  =  A  +Bx  +Cx*  +  Dx 3 

where  A,  B,  C,  and  D  are  specified  constants.  A  similar  temperature  distribution  with  a  modified 
value  for  A  was  imposed  near  the  centerline  of  the  rod  at  t  =  0  so  that,  at  steady  state,  the  wall  tem¬ 
perature  distribution  agreed  with  the  experimental  curve  fit.  The  steady-state  temperature  field  in  the 
rod  and  in  the  bulk  fuel  was  then  determined  and  stored.  Since  the  results  were  similar  to  the  solution 
shown  above,  they  are  not  presented  here.  Subsequently,  the  calibration  studies  were  carried  out.  In 
the  following,  we  describe  some  of  the  more  pertinent  results  of  the  calibration  studies. 

To  determine  the  form  of  the  deposit  growth  equation  (68),  quite  a  few  parameters  had  to  be 
determined.  For  the  first  mechanism,  i.e.,  the  catalytic  wall  reaction,  two  parameters  had  to  be  deter¬ 
mined,  Aa  and  Ew.  From  the  JFTOT  data  (UTC,  1988),  it  had  appeared  that  the  wall  activation 
energy  had  a  value  in  the  range  of  35-40  kcal/mole  for  the  types  of  JP-5  fuel  studied.  The  experimen¬ 
tal  data  for  the  nonlinear  growth  of  the  deposit  was  for  the  JP-5  fuel  with  an  activation  energy  of 
Ew  =  39  kcal/mole.  We  therefore  decided  to  fix  the  wall  activation  energy  to  this  value.  Thus,  the 
calibration  of  the  linear  growth  of  the  deposit  due  to  the  catalytic  wall  reaction  required  the  determi¬ 
nation  of  the  collision  frequency,  i.e.,  Av 

For  the  second  term,  i.e.,  the  nonlinear  growth  of  the  wall  deposit  due  to  the  random  impaction 
and  adhesion  of  the  solid  particles  formed  due  to  phase  change  of  the  liquid  precursor  near  the  wall, 
we  had  to  determine  the  activation  energies,  Ev  and  EL,  and  the  collision  frequencies,  APand  AL,  for 
the  bulk  fuel  precursor  production  and  loss  mechanisms,  Equations  (59)  and  (60).  We  followed  the 
study  of  Krazinski  et  al.  (1989)  and  used  EP  =  30-40  kcal/mole  and  EL  =  35-45  kcal/mole  for  the 
model.  Most  of  the  results  shown  here  employed  EP-  30  kcal/mole  and  EL  =  35  kcal/mole.  The  con¬ 
stant  Bw  in  Equation  (68)  was  set  to  unity  for  lack  of  any  additional  experimental  information  on  this 
mechanism.  Then  the  calibration  process  required  the  determination  of  the  constants  Ap  and  AL.  The 
details  of  all  the  calibration  calculations  carried  out  will  not  be  discussed  here;  instead,  we  will  focus 
on  some  characteristic  results. 

Figure  5  illustrates  the  time-dependent  growth  of  the  solid  deposit  on  the  wall  for  a  variety  of 
choices  of  the  calibration  constants.  All  curves  are  identified  in  the  figure.  With  Bw  =  0  and  C„  =  0, 
we  obtain  the  linear  growth  of  the  deposit  due  to  the  wall  catalytic  reaction  with  any  value  of  Av 
The  calibration  study  indicated  that  Am  =  104  m/s  is  a  good  choice.  This  figure  also  shows  the  effect 
of  just  the  nonlinear  second  term  with  Aw  =  I0-8  m/s  and  C„  =  0  by  another  curve.  The  nonlinear 
term  allows  the  deposit  to  grow  very  slowly  at  first  and  then,  when  the  liquid  precursor  begins  to  form 
in  the  bulk  fuel  next  to  the  wall,  the  nonlinear  term  rapidly  grows  in  importance.  Finally,  this  figure 
also  shows  the  effect  of  combining  the  various  mechanisms  for  the  deposit  growth.  Comparison  with 
the  experimental  data  shows  that  values  of  AP=  9xl07sec-1  and  AL  =  9xl08sec-1  for  the  chosen 
values  of£>  and  EL  appear  to  provide  close  agreement  with  the  experimental  data. 
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Figure  5. 


Time-Dependent  Growth  of  the  Deposit  at  Peak  Temperature  Location  for 
Various  Models.  Jet  fuel  is  NACP  FF-8  (JP-5);  peak  experimental  temperature  is 
550®F  (286°C).  Model  activation  energies:  Ew  ■  39  kcal/mole,  Ep  =  30  kcal/mole, 
£l  -  35  kcal/mole.  Also,  unless  otherwise  noted  in  the  figure,  Aw  ■  104  m/s, 

Ap  »  9xl07  sec*1.  At  «  9x10®  sec'1,  Bm  -  1. 
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Figure  6a  shows  the  axial  variation  of  the  solid  deposit  distribution  for  various  times  for  a  particu¬ 
lar  calibration  study.  The  rapid  nonlinear  growth  of  the  deposit  with  increase  in  time  can  be  seen  in 
this  figure.  Figure  6b  shows  solid  deposit  distribution  at  the  last  time  on  which  the  wall  temperature 
is  also  plotted.  The  deposit  thickness  increases  with  increasing  temperature,  with  the  maximum  depo¬ 
sition  occurring  near  the  peak  temperature  on  the  wall. 

Figure  7a  shows  the  growth  of  the  liquid  precursor  (by  mass  fraction  Y)  near  the  edge  of  the  solid 
deposit  (at  r  =  r2)  for  a  representative  case.  It  can  be  seen  that  peak  production  of  the  liquid  deposit 
precursor  also  occurs  in  the  region  where  the  temperature  is  high.  As  time  increases,  however,  the 
deposit  precursor  distribution  near  the  rod  surface  becomes  more  uniform  in  the  axial  direction.  In 
the  radial  direction,  the  temperature  drops  rapidly  and  the  net  liquid  deposit  precursor  production 
also  drops  rapidly. 

Figure  7b  shows  the  axial  variation  of  the  temperature  of  the  bulk  fuel  adjacent  to  the  solid  depo¬ 
sit  during  the  calculation  for  a  case  with  the  first  two  mechanisms  of  deposit  growth  [Equation  (68)]. 
As  noted  in  the  formulation,  we  allowed  the  temperature  of  the  wall  and  the  bulk  fuel  to  adjust  to  the 
growth  effects  of  the  deposit  by  virtue  of  the  boundary  matching  described  in  Section  3.2.  Since  we 
had  assumed  the  thermal  conductivity  of  the  deposit  to  be  five  times  smaller  than  that  of  the  fuel  (i.e., 
Ad  =  0.2),  this  boundary  matching  results  in  an  increase  in  the  wall  temperature  with  time  as  the 
deposit  grows  and  a  corresponding  decrease  in  the  bulk  fuel  temperature  with  growth  in  the  deposit. 
This  is  due  to  the  inhibiting  effect  of  the  growing  deposit  layer  on  the  heat  transfer  process.  Figure  7b 
shows  that  as  the  deposit  grows,  the  temperature  of  the  bulk  fuel  starts  to  decrease,  especially  in  the 
regions  where  the  deposit  is  maximum.  As  modeled  here,  this  phenomena  appears  to  agree  with  the 
observations  of  Marteney  (1988).  The  wall  temperature  also  begins  to  increase;  however,  the  effect  of 
the  solid  deposit  on  the  temperature  is  not  very  significant  in  the  range  of  growth  modeled  here. 
Therefore,  as  a  first  approximation,  we  could  have  neglected  this  effect  and  assumed  that  the  deposit 
growth  does  not  affect  the  heat  transfer  between  the  heated  surface  and  the  fuel.  Additional  study 
using  different  values  for  the  thermal  conductivity  (for  example,  Ad  =  0.02)  showed  that  in  this  case 
the  effect  of  deposit  growth  on  the  heat  transfer  process  could  be  substantial.  The  significance  of  this 
mechanism  can  only  be  quantified  if  the  physical  properties,  such  as  the  dependence  of  the  deposit 
thermal  conductivity  on  temperature,  kj(T),  are  known. 

This  validation  study  did,  however,  allow  us  to  calibrate  the  model  for  the  nonlinear  variation 
observed  in  the  JFTOT.  Subsequent  to  this  study  (near  the  end  of  Phase  I),  we  obtained  the  data  of 
Marteney  and  Spadaccini  (1986)  and  the  computed  calibration  by  Krazinski  et  al.  (1990).  This  data 
was  obtained  using  a  different  test  apparatus  and  included  data  on  deposit  growth  which  showed  that 
the  wall  activation  energy  was  probably  much  lower  than  that  estimated  from  the  JFTOT  data.  A 
choice  of  £.  =  8  kcal/mole  was  considered  more  approoriate.  Since  this  wais  quite  different  from  the 
earlier  choice  of  £„  =  39  kcal/mole  used  in  our  calibration  study  described  above,  we  decided  to 
study  the  effect  of  changing  the  wall  activation  energy. 

Note  that  our  model  had  two  primary  mechanisms,  the  wall  catalytic  reaction  and  the  nonlinear 
growth  due  to  the  bulk  fuel  precursor  undergoing  phase  transition  and  impacting  in  a  random  sense  on 
the  wall.  Thus,  if  our  model  was  formulated  correctly,  and  if  the  wall  activation  energy  £tt  was 
reduced  to  8  kcal/mole,  then  the  only  calibration  that  should  be  required  is  to  determine  the  new 
value  of  Aa  with  Bw  =  0  and  C„  =  0.  Then  if  the  new  calibrated  wall  catalytic  reaction  mechanism  is 
combined  with  the  model  for  the  bulk  fuel  phase  change  calibrated  earlier,  the  predicted  growth  of 
the  deposit  with  time  should  still  agree  with  the  experimentally  observed  variation.  This  was  in  fact 
observed  in  our  study  as  shown  below. 
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a.  As  a  function  of  time 


«v 
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b.  After  approximately  4.5  hours  (numbers  on  curve  indicate  wail  temperature) 

Figure  6.  Axial  Variation  of  the  Deposit  Thickness  on  the  Rod  Surface.  Test  conditions: 
Ew  -  39  kcal/mole,  Ep  -  30  kcal/raole,  El  *  35  kcal/mole,  Aw  -  104  m/s, 

Ad  -  9xl07  sec-1,  Al  -  9x10*  sec-1,  Cw  -  1,  Bw  -  1. 
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Figure  8  shows  the  time  growth  of  the  deposit  on  the  wall  using  the  new  calibration  for  the  wall 
activation  energy  of  Ew  =  8  kcal/mole.  As  expected,  once  the  linear  growth  model  for  the  wall  cata¬ 
lytic  conversion  was  modeled  (i.e.,  Aw  was  determined),  the  combined  model  again  predicts  the  same 
nonlinear  growth.  This  appears  to  indicate  that  the  proposed  nonlinear  mechanism  may  have  some 
real  validity.  Of  course,  the  validity  of  the  model  can  only  be  confirmed  if  additional  experimental 
data  is  obtained  to  determine  the  proposed  Vh  variation  for  the  growth  of  the  deposit  thickness  with 
time  does  in  fact  occur  in  the  actual  experiments.  Future  experiments  could  be  designed  to  determine 
if  the  mechanism  proposed  in  this  model  is  correct 

Figures  9a  and  9b  show  the  axial  growth  of  the  deposit  with  time  for  the  cases  in  which  the  third 
mechanism  was  neglected  and  for  the  case  in  which  the  modified  form  of  the  third  mechanism  [Equa¬ 
tion  (81)]  was  included,  respectively.  Comparison  shows  that  for  the  case  in  which  the  third  mechan¬ 
ism  was  neglected,  the  axial  distribution  is  quite  smooth.  However,  with  the  modified  mechanism,  the 
axial  growth  of  the  deposit  downstream  of  the  peak  temperature  location  is  higher.  This  trend  appears 
to  be  similar  to  the  increase  in  the  deposit  thickness  downstream  cr  th_  peak  temperature  observed  in 
the  JFTOT  data  (UTC,  1988).  In  Figure  9c,  the  temperature  is  plotted  on  the  curve  for  the  axial  vari¬ 
ation  of  the  deposit  growth  at  the  end  of  the  calculation  in  which  the  third  mechanism  was  included. 
For  the  same  wall  temperature,  the  thickness  of  the  deposit  downstream  of  the  peak  temperature  is 
larger  when  compared  to  the  upstream  value. 


6.  CONCLUSIONS  AND  RECOMMENDATIONS 

In  this  Phase  I  study,  a  detailed  mathematical  model  was  developed  that  would  allow  the  modeling 
of  the  various  possible  physical  mechanisms  that  lead  to  the  thermal  degradation  of  jet  fuels  and  the 
subsequent  formation  of  solid  deposits  on  a  heated  surface.  The  mathematical  formulation  would 
allow  the  calculation  of  the  heat  conduction  through  the  heated  metal  surface,  the  heat  transfer 
through  the  solid  deposit,  the  conservation  of  mass  momentum  and  energy  in  the  bulk  fuel,  the  chemi¬ 
cal  kinetic  processes  in  the  bulk  fuel,  and  the  time-dependent  mechanism  for  the  growth  of  the  solid 
deposit  on  the  wall.  In  general,  the  model  developed  is  applicable  to  three-dimensional  flows  of  fuels 
in  appropriate  fuel  systems.  If  a  turbulence  model  is  included,  then  turbulent  flows  can  also  be  stu¬ 
died.  Time-  and  space-accurate  numerical  schemes,  both  explicit  and  implicit,  were  developed  for  use 
with  this  model  in  this  study. 

Although  the  mathematical  model  formulated  was  quite  general,  it  was  not  possible  to  solve  the 
complex  mathematical  equations  for  a  multispecies  two-phase  flow  in  Phase  I.  Therefore,  a  simplified 
model  was  reduced  from  the  general  formulation  using  axisymmetric  assumptions  and  then  solved 
using  the  numerical  scheme  developed  in  this  study.  The  steady-state  temperature  distribution  in  the 
heating  dement  and  in  the  bulk  fuel  was  calculated  using  this  numerical  scheme.  It  was  demonstrated 
that  under  the  axisymmetric  approximation,  the  temperature  distribution  will  not  agree  with  the 
observed  JFTOT  variation  if  the  constant  heat  flux  boundary  condition  is  used.  By  modifying  the 
heat  flux  distribution,  however,  we  were  able  to  reproduce  the  JFTOT  variation.  Since  this  was  as 
good  as  using  the  curve  fits  for  the  experimental  data,  we  used  the  curve  fits  to  obtain  the  steady-state 
temperature  field  in  the  bulk  fuel  for  the  calibration  studies. 

With  the  computed  temperature  field,  we  carried  out  a  series  of  calibration  studies  to  obtain  the 
calibration  constants  for  the  deposit  growth  and  the  liquid  precursor.  These  studies  showed  that  the 
proposed  growth  mechanism  could  be  made  to  agree  with  the  JFTOT  experimental  data  on  the  non¬ 
linear  growth  rate  (UTC,  1988)  quite  well.  This  appears  to  indicate  that  the  proposed  mechanism  may 
have  some  validity. 
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Figure  8.  Time-Dependent  Growth  of  Deposit  Thickness  at  the  Peak  Temperature  Location 
for  Ew  *  8  Kcal/mole.  Jet  fuel  and  peak  experimental  temperature  same  as  in 
Fig.  5.  Also,  unless  otherwise  noted  in  the  figure,  £p  «  30  kcal/mole,  £l  *b  35 
kcal/mole,  Am  ■  8xl0‘9  m/s,  Ap  «  OxlO^ec*1,  A l  -  OxlO^ec*1,  Bm  ■  1,  Cw  =  0. 


a.  As  a  function  of  time.  Test  conditions:  £«■!  kcal/mole, 
£p  -  30  kcal/mole,  £l  -  35  kcal/mole,  Aw  -  8x10'®  m/s, 

A p  -  9xl07  sec*1,  A l  -  9x10*  sec*1,  Bw  -  1,  Cw  -  1. 


Figure  9.  Axial  Variation  of  Deposit  Thickness 
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b.  As  a  function  of  time.  Test  conditions  same  as  in  Fig.  9a  except  Cw  »  2.5. 


CM 


c.  After  approximately  4.5  hours  (numbers  on  curve  indicate  wall  temperature). 
Test  conditions  same  as  in  Fig.  9b. 

Figure  9.  Axial  Variation  of  Deposit  Thickness  (Cont.) 
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One  of  the  major  hurdles  facing  the  validation  of  CFD  models  such  as  the  one  developed  here  to 
study  the  thermal  degradation  process  is  that  many  of  the  mechanisms  leading  to  the  growth  of  the 
solid  deposit  are  largely  unknown.  In  addition,  the  experimental  data  available  are  quite  limited  both 
in  quality  and  quantity.  In  the  present  study,  we  used  primarily  one  set  of  data  based  on  the  JFTOT 
experiments  to  calibrate  our  model.  The  reduced  model  used  in  Phase  I  appears  to  show  good  agree¬ 
ment  with  the  experimental  data  of  the  JFTOT  study  (UTC,  1988).  However,  it  is  not  clear  if  the 
JFTOT  data  itself  is  an  accurate  representation  of  the  actual  thermal  degradation  process.  In  addi¬ 
tion,  the  data  makes  no  attempt  to  quantify  the  physical  process  leading  to  deposit  growth  and  pro¬ 
vides  only  gross  information  on  the  amount  of  deposit  on  the  heated  surface.  Thus,  although  the 
model  appears  to  be  validated  by  comparison  with  the  data,  it  is  not  clear  if  the  physical  mechanism 
proposed  for  the  growth  of  the  deposit  is  a  valid  representation  of  the  actual  physics  of  the  complex 
chemical  processes  involved  in  the  thermal  degradation  of  the  jet  fuels.  Therefore,  it  is  probably 
premature  to  further  extend  this  model  at  this  stage  to  include  the  two-phase  full  model  formulated  in 
this  study.  More  controlled  and  detailed  experimental  studies  must  be  conducted  before  additional 
evaluation  of  CFD  models  can  be  carried  out.  Such  experiments  can  be  considered  in  the  next  phase 
so  that  additional  data  is  available  before  further  model  development  is  carried  out. 
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